//VST\ 


MIT/WHOI  97-23 


Massachusetts  Institute  of  Technology 
Woods  Hole  Oceanographic  Institution 


Joint  Program 
in  Oceanography/ 
Applied  Ocean  Science 
and  Engineering 


1930 


DOCTORAL  DISSERTATION 


Modeling  of  Monostatic  Bottom  Backscattering  from 
Three-Dimensional  Volume  Inhomogeneities  and 
Comparisons  with  Experimental  Data 

by 


iDlSTHibUTION  STA' 


Approved  for  public  releosi^ 
Distribution  Unlimited 


Dan  Li 


eptember  1997 


19980505  058 


0 


% 


MITAVHOI 

97-23 


Modeling  of  Monostatic  Bottom  Backscattering  from  Three-Dimensional  Volume 
Inhomogeneities  and  Comparisons  with  Experimental  Data 


by 

Dan  Li 


Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts  02139 

and 

Woods  Hole  Oceanographic  Institution 
Woods  Hole,  Massachusetts  02543 

September  1997 

DOCTORAL  DISSERTATION 


Funding  was  provided  by  Office  of  Naval  Research  Contract  N0001 4-91 -J-1 379. 


Reproduction  in  whole  or  in  part  is  permitted  for  any  purpose  of  the  United  States  Government.  This 
thesis  should  be  cited  as:  Dan  Li,  1997.  Modeling  of  Monostatic  Bottom  Backscattering  from 
Three-Dimensional  Volume  Inhomogeneities  and  Comparisons  with  Experimental  Data  Ph  D  Thesis 

MIT/WHOI,  97-23. 


Approved  for  publication;  distribution  unlimited. 


Approved  for  Distribution: 


Timothy  K.  Stanton,  Chair 

Department  of  Applied  Ocean  Physics  and  Engineering 


Paola  Malanoppe-Rizzoli 

MIT  Director  of  Joint  Program 


WHOI  Dean  of  Graduate  Studies 


Modeling  of  Monostatic  Bottom  Backscattering  from 
Three-Dimensional  Volume  Inhomogeneities  and 
Comparisons  with  Experimental  Data 

by 

Dan  Li 

B.  S.,  Electrical  Engineering 
University  of  Science  and  Technology  of  China  (1992) 

M.  S.,  O.  E.,  Ocean  Engineering 
M.I.T.  /  Woods  Hole  Oceanographic  Institution  (1995) 

Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 

DOCTOR  OF  PHILOSOPHY  IN  OCEANOGRAPHIC  ENGINEERING 

at  the 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 

and  the 

WOODS  HOLE  OCEANOGRAPHIC  INSTITUTION 
September  1997 

©  Dan  Li,  1997.  All  rights  reserved. 

The  author  hereby  grants  to  MIT  and  WHOI  permission  to  reproduce  and 
to  distribute  copies  of  this  thesis  document  in  whole  or  in  part. 


Signature  of  Author 


Certified  by 


Certified  by 


Accepted  by 


Department  of  Ocean  Engineering,  MIT  and  the 
MIT/WHOI  Joint  Program  in  Oceanographic  Engineering 


raphic  Engineering 


y  Dr.  George  V.  Frisk 
Senior  Scientist,  Woods  Hole  Oceanographic  Institution 


Associate  Scientist,  Woods  Ho^ 


^r^hic  Institution 
Thesis  Supervisor 


Chairman,  Joint  Committee  for 


'^Pro^^ofiHenrik  Schmidt 
l&eatlographic  Engineering, 


Massachusetts  Institute  of  Technology- Woods  Hole  Oceanographic  Institution 


1 


Modeling  of  Monostatic  Bottom  Backscattering  from 
Three-Dimensional  Volume  Inhomogeneities  and 
Comparisons  with  Experimental  Data 

by 

Dan  Li 


Submitted  to  the  Massachusetts  Institute  of  Technology/  Woods  Hole 
Oceanographic  Institution  in  partial  fulfillment  of  the  requirements  for  the  degree  of 

DOCTOR  OF  PHILOSOPHY  IN  OCEANOGRAPHIC  ENGINEERING 


Abstract 

Acoustic  propagation  in  the  ocean  inevitably  encounters  inhomogeneities  of  various 
types,  which  give  rise  to  scattering.  Acoustic  scattering  from  rough  water/bottom 
interfaces  comprised  of  exposed  rocks  and  sea  mountains  gives  way  to  volumetric 
scattering  in  areas  with  flat  interfaces  and  thick  sediment  cover.  The  data  analysis 
of  the  ARSRP  backscattering  experiment  revealed  that  random  inhomogeneities  in 
two  irregular  layers  beneath  the  seafloor  were  the  primary  contributors  to  oblique 
backscattering  in  a  sediment  pond  on  the  western  flank  of  the  Mid- Atlantic  Ridge.  In 
this  thesis,  an  attempt  has  been  made  to  model  monostatic  backscattering  from  3-D 
volume  inhomogeneities  in  the  sediment  and  to  compare  the  results  with  the  ARSRP 
backscattering  data. 

A  scattering  process  cannot  be  modeled  correctly  without  a  proper  account  of  the 
incident  field.  Several  approximate  propagation  models  have  been  evaluated  against 
the  exact  solution,  while  the  appropriateness  of  using  the  equivalent  surface  scatter¬ 
ing  strength  in  volume  scattering  characterizations  is  studied.  This  study  concludes 
that  precautions  need  to  be  taken  in  modeling  both  the  propagation  effects  and  the 
scattering  mechanisms  associated  with  the  bottom  volume  scattering  process. 

A  volume  scattering  model  based  on  perturbation  theory  and  the  Born  approxi¬ 
mation  is  developed  incorporating  contributions  from  both  sound  speed  and  density 
fluctuations.  With  the  propagation  part  handled  accurately  by  OASES  and  random 
fluctuations  generated  effectively  by  a  new  scheme  modified  from  the  spectral  method, 
the  model  is  capable  of  simulating  the  monostatic  backscattered  field  and  time  series 
due  to  3-D  volumetric  sediment  inhomogeneities.  Both  the  characteristic  length  scale 
and  power  spectrum  descriptions  of  the  random  inhomogeneities  are  shown  to  have 
great  impact  on  the  backscattered  field  by  parameter  studies  in  a  free-space  scenario. 
The  important  roles  played  by  horizontal  anisotropy  and  the  vertical  correlation  of 
the  random  field  have  been  demonstrated.  Density  fluctuations  are  further  confirmed 
to  be  the  dominant  force  in  backscattering.  The  model  matches  the  ARSRP  backscat¬ 
tering  data  very  well,  with  the  fluctuations  of  sound  speed  and  density  in  the  two 
irregular  layers  described  by  a  power  law  type  of  power  spectrum. 
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Chapter  1 


Introduction 


1.1  General  Background 

Acoustic  propagation  in  the  ocean  will  inevitably  encounter  various  types  of  inho¬ 
mogeneities.  Rough  sea  surfaces  and  seafloors,  air  bubbles,  schools  of  fish,  internal 
waves,  and  variations  of  sediment  properties  constitute  only  a  small  part  of  a  wide 
spectrum.  All  of  these  inhomogeneities  have  the  capability  of  reradiating  a  certain 
amount  of  incident  acoustic  energy,  the  process  of  which  is  called  scattering. 

With  abilities  to  affect  the  distribution  of  the  acoustic  energy,  scattering  will  un¬ 
doubtedly  play  a  significant  role  in  sonar  operation,  which  remains  a  primary  means 
to  explore  the  ocean.  Analogous  to  optical  scattering,  which  makes  the  world  vivid  to 
us,  acoustic  scattering  enables  us  to  “see”  the  underwater  world.  Broad  applications 
include  target  detection  and  localization,  geological  surveys  of  the  ocean  seafloor,  and 
monitoring  of  environmental  changes,  to  name  a  few.  On  the  other  hand,  contrary  to 
its  above  role  as  signal,  the  reverberation  can  act  as  background  noise  to  operations 
such  as  acoustic  navigation  and  communication.  Especially  since  it  is  induced  by 
the  transmitted  signal  itself,  it  would  be  extremely  difficult  to  achieve  a  successful 
elimination  without  the  knowledge  of  the  excitation  mechanisms.  Needless  to  say, 
understanding  scattering  mechanisms  is  essential  to  sonar  operations  regardless  of 
the  role  scattering  is  playing. 
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Usually,  acoustic  scattering  is  categorized  into  rough  surface  scattering  and  vol¬ 
ume  scattering.  For  many  years,  scattering  from  rough  surfaces  has  been  the  center  of 
attention.  Rough  sea  surfaces,  exposed  rocks  on  the  seafloor,  and  large  topographical 
features  such  as  sea  mountains  and  ridges  are  considered  the  dominant  contributors  to 
bottom  scattering.  On  the  other  hand,  bottom  volume  scattering  has  attracted  more 
and  more  interest  in  recent  years.  Although  less  obvious  and  less  well  understood, 
more  evidence  shows  that  the  inhomogeneities  within  the  sediment  can  be  a  dominant 
factor  in  total  bottom  scattering,  especially  when  the  seafloor  is  relatively  flat  and/or 
the  the  bottom  attenuation  is  small.  In  his  early  work  on  the  measurement  of  the 
bottom  backscattering  strength,  Merklinger  [2]  speculated  that  subbottom  inhomo¬ 
geneities  may  be  primary  contributors  to  the  backscattered  field.  Jackson  et  al.  [3,  4] 
and  Lyons  et  al.  [5]  also  recognized  the  importance  of  volumetric  inhomogeneities 
in  the  bottom  scattering.  Tang  et  al.’s  [6]  analysis  of  high-frequency  scattering  data 
showed  that  the  gas  voids  in  the  sediment  are  probable  significant  scatterers. 

In  order  to  understand  a  physical  process  such  as  scattering,  we  will  have  to  re¬ 
sort  to  both  experiment  and  modeling.  Due  to  various  reasons,  sediment  volume 
scattering  is  difficult  to  study.  First  of  all,  it  is  extremely  difficult  to  directly  mea¬ 
sure  the  properties  of  the  sediment.  As  a  result,  a  lot  of  modeling  work  is  based  on 
very  limited  knowledge  of  the  ground  truth.  Second,  it  is  no  easy  task  designing  a 
good  scattering  experiment  as  well.  In  the  usual  deep-water  scattering  experiments, 
the  source  and  receivers  are  close  to  the  sea  surface.  Consequently,  the  insonified 
region  is  so  large  that  the  strong  returns  from  sea  mountains  are  likely  to  mask  the 
weak  returns  from  the  sediment  volume.  In  shallow  water,  the  waveguide  effects  will 
inevitably  complicate  the  scattered  field,  making  it  arduous  to  pinpoint  the  volume 
scattering  effect  [7,  8].  Third,  the  conventional  single  source/single  receiver  configu¬ 
ration  is  unable  to  distinguish  the  volumetric  scattered  returns  in  different  directions, 
which  may  result  in  misinterpretation  of  the  data.  To  avoid  all  of  the  above  problems, 
3-n  ideal  experimental  scenario  would  be  to  have  a  deep-towed  acoustic  source  and 


18 


receiving  array  collecting  data  over  a  deep-sea  sediment  pond  area,  where  the  acous¬ 
tic  propagation  is  relatively  simple  in  the  absence  of  the  waveguide  effects  manifested 
in  shallow  water  and  the  scattered  returns  are  primarily  from  the  volume  inhomo¬ 
geneities  in  the  sediment  because  the  deep-sea  water/sediment  interfaces  are  normally 
flat.  Fortunately  for  us,  the  experiment  to  be  described  next  fit  this  scenario  very  well. 


1.2  Data  Analysis  Results  of  ARSRP  Backscatter- 
ing  Experiment 

As  part  of  the  Office  of  Naval  Research  Bottom/Subbottom  Acoustic  Reverberation 
Special  Research  Program  (ARSRP),  a  backscattering  experiment  was  conducted  over 
a  sediment  pond  on  the  western  flank  of  the  Mid- Atlantic  Ridge  in  July  1993.  A  chirp 
source  and  vertical  line  array  were  deep-towed  200  m  to  400  m  above  the  seafloor, 
while  the  sediment  thickness  was  up  to  400  m  in  the  middle  of  the  sediment  pond. 
For  a  transmitted  acoustic  source  level  of  approximately  200  dB  re  l^uPa  @  Im  and 
a  source  frequency  of  250-650  Hz,  the  acoustic  signal  can  penetrate  to  the  rock  base¬ 
ment  beneath  the  sediment  layer  in  spite  of  sediment  attenuation  and  geometrical 
spreading.  Details  of  the  experiment  description  are  presented  in  Chapter  6.  The 
following  are  the  data  processing  results  from  Ref.  [9,  10]. 

Taking  advantage  of  the  vertical  line  array,  endfire  beamforming  yields  signals 
which  are  dominated  by  normal  incidence  returns.  With  the  ship  moving  across  the 
sediment  pond  area,  sediment  profiling  can  be  obtained  by  aligning  the  normal  in¬ 
cidence  returns  with  respect  to  reflection  arrivals  at  the  water/sediment  interface. 
After  employing  an  edge  detection  algorithm  [11],  the  inferred  sediment  structure  is 
shown  in  Fig.  1-1  at  the  east  side  of  the  sediment  pond.  Evidently,  layering  with 
gentle  horizontal  changes  is  the  main  characteristic  of  the  sediment.  However,  some 
random  features  can  be  seen  in  two  irregular  layers  beneath  the  water/sediment  in- 
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Figure  1-1:  The  inferred  sediment  structure  at  the  east  side  of  the  sediment  pond. 


terface,  from  0.02s  to  0.045s  and  from  0.065s  to  0.11s.  A  geological  interpretation 
suggests  that  they  are  most  likely  to  be  turbidity  layers  [12].  What  interests  us  is  to 
find  out  whether  the  inhomogeneities  contribute  to  backscattering  in  oblique  direc¬ 
tions.  After  subtracting  the  returns  in  normal  and  near  normal  incidence  directions, 
the  backscattered  returns  at  oblique  grazing  angles  can  then  be  isolated  by  multi¬ 
ple  constraints  beamforming.  We  average  the  envelope  of  the  beamformed  results 
over  8  consecutive  pings  and  further  smooth  the  curve  by  a  low-pass  filter  in  order 
to  detect  any  potential  trend.  The  averaged  envelope  is  shown  in  Fig.  1-2.  Three 
groups  of  peaks  appear  in  each  look  direction.  Assuming  that  these  peaks  occur  at 
horizontal  interfaces  and  using  a  nominal  sound  speed  of  1530m/s,  we  have  the  best 
fit  indicated  by  the  solid  lines.  Here  time  0  corresponds  to  the  arrival  time  of  the 
first  water /sediment  interface  refiection.  From  the  travel  times,  we  find  that  the  first 
two  imaginary  interfaces  would  reside  at  the  locations  of  the  above-mentioned  two 
irregular  layers,  where  actual  interfaces  don’t  appear  to  exist.  We  gather  that  the  ob¬ 
served  strong  backscattering  is  probably  caused  by  volumetric  inhomogeneities.  The 
third  group  is  associated  with  sediment/basement  interface  according  to  travel  time. 
Therefore,  we  tentatively  conclude  that  the  backscattered  field  in  oblique  directions 
is  due  to  sediment  inhomogeneities  in  the  two  irregular  layers. 

The  same  procedure  is  executed  for  pings  collected  at  the  west  side  of  the  sed¬ 
iment  pond.  Figure  1-3  shows  the  sediment  structure.  Interestingly,  the  upper  ir¬ 
regular  layer  is  absent  on  this  side  of  the  sediment  pond.  Therefore,  if  the  above 
hypothesis  that  inhomogeneities  in  the  irregular  layers  are  the  main  contributors  to 
the  oblique  backscattering  is  true,  one  would  see  only  one  group  of  peaks  correspond¬ 
ing  to  the  lower  irregular  layer  instead  of  two  groups  of  peaks  at  the  beginning  of 
the  beamformed  output  record.  Figure  1-4  shows  the  average  envelope  of  oblique 
angle  beamforming  results.  The  striking  point  here  is  that  the  first  group  of  peaks 
does  disappear,  which  confirms  our  finding  that  the  irregular  layers  beneath  the  wa¬ 
ter/sediment  interface  are  the  dominant  contributors  to  backscattering  in  oblique 
directions.  Notice  that  there  is  not  a  group  of  peaks  close  to  time  0.02s  at  60  degree 
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Figure  1-3:  The  inferred  sediment  structure  at  the  west  side  of  the  sediment  pond. 
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grazing  angle  in  both  Fig.  1-2  and  Fig.  1-4,  which  suggests  the  insignificance  of 
backscattering  at  the  water/sediment  interface.  This  has  been  independently  con¬ 
firmed  by  Jeffrey  and  Ogden  [13]. 


1.3  Objective  and  Approach 

The  results  of  experimental  data  analysis  are  encouraging.  Now  it  is  up  to  the  model¬ 
ing  to  help  achieve  better  understanding  of  the  scattering  mechanisms.  As  a  matter  of 
fact,  the  objective  of  this  thesis  is  to  develop  a  monostatic  backscattering  model  that 
can  incorporate  3-D  bottom  volumetric  inhomogeneities  and  to  compare  the  results 
with  the  ARSRP  backscattering  data. 

Since  the  propagation  of  sound  between  the  source/receiver  and  the  scatterer 
should  be  treated  as  an  integral  part  of  the  scattering  modeling,  as  the  first  step, 
efforts  will  be  made  to  evaluate  propagation  approximations  incorporated  in  Hines’ 
[14],  Ivakin’s  [15]  and  Mourad  and  Jackson’s  [1]  scattering  models.  Comparisons 
with  the  full-wave  solution  obtained  by  a  wavenumber  integration  method  [16]  are 
expected  to  show  the  validity  and  limitations  of  those  approaches.  The  effectiveness 
of  the  equivalent  surface  scattering  strength  in  characterizing  the  volume  scattering 
process  is  going  to  be  examined  as  well. 

Due  to  the  fact  that  the  equivalent  surface  scattering  strength  is  inappropriate  in 
describing  the  angular  dependence  of  the  volume  scattering  when  multi-path  effects 
cannot  be  ignored,  we  want  the  model  to  possess  an  ability  to  overcome  this.  The 
most  direct  way  to  conduct  a  model/data  comparison  of  bottom  volume  scattering 
is  to  compare  them  in  the  time  domain.  Scattering  from  volume  inhomogeneities  is 
usually  considered  as  a  random  process,  it  is  therefore  practically  impossible  to  syn¬ 
thesize  time  series  which  can  match  a  single  trace  of  the  scattered  field.  Instead,  one 
would  try  to  compare  the  average  shape,  or  frequency  components,  of  the  synthesized 
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time  series  with  that  of  the  experimental  data.  A  forward  model  that  can  incorpo¬ 
rate  random  inhomogeneities  and  generate  reverberation  time  series  would  definitely 
improve  our  capability  of  simulating  scattering  experiments. 

In  the  model,  the  calculation  of  the  incident  field  will  be  handled  by  a  numerical 
wavenumber  integration  technique  using  OASES  [16].  This  would  provide  the  model 
with  versatility  to  deal  with  scattering  involving  complicated  propagation  scenario 
such  as  a  sound  speed  gradient  and  even  waveguide  effects  in  shallow  water.  Similar 
to  Morse  and  Ingard’s  [17]  and  Tang’s  [8]  approach,  a  perturbation  method  will  be 
applied  to  determine  the  scattered  part  of  the  field.  Here,  sound  speed  and  den¬ 
sity  fluctuations  in  the  fluid  sediment  are  considered  the  cause  of  scattering.  Rough 
surface  scattering  will  not  be  included  here  but  can  be  added  in  a  straightforward 
manner.  While  the  necessity  to  include  poro-elastic  effects  in  a  sediment  scattering 
model  remains  an  open  question  [18],  a  fluid  model  is  appropriate  for  the  soft  sed¬ 
iments  and  low  frequencies  under  consideration.  The  volume  inhomogeneities,  i.e., 
the  sound  speed  and  density  variations,  are  generated  by  the  spectral  method  [19]. 
The  random  field  is  assumed  to  be  correlated  in  three  dimensions  and  defined  by  a 
power  spectrum.  By  taking  advantage  of  the  fact  that  the  source  and  receiver  are 
both  omnidirectional  and  co-located  in  monostatic  backscattering  experiment  scenar¬ 
ios,  the  usually  computationally  intensive  procedure  can  be  carried  out  effectively. 
The  model  is  expected  to  interpret  the  phenomena  which  have  been  observed  in  the 
ARSRP  site  A  sediment  scattering  experiment  [9,  10]. 


1.4  Literature  Review 

The  modeling  of  scattering  from  both  rough  interfaces  and  volume  inhomogeneities 
has  been  an  active  research  area  for  many  years,  and  the  related  literature  is  vast. 
We  will  present  a  brief  overview  of  it. 
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A  large  amount  of  work  has  been  reported  studying  acoustic  scattering  from  rough 
surfaces.  There  are  two  well-known  theoretical  methods  for  calculating  acoustic  scat¬ 
tering  from  rough  surfaces:  the  Rayleigh-Rice  method  and  the  Kirchhoff  method. 
Based  on  the  small  perturbation  approximation,  the  Rayleigh-Rice  method  treats 
the  rough  surface  as  a  perturbation  to  a  smooth  plane  and  the  resulting  scattering 
coefficient  due  to  the  presence  of  roughness  is  calculated.  The  Kirchhoff  approxima¬ 
tion  approach  is  more  pertinent  to  large-scale  roughness,  and  the  scattering  surface 
is  assumed  to  be  sufficiently  smooth  so  that  the  tangent  plane  at  any  point  of  the 
surface  determines  the  reflection  properties.  Some  efforts  are  also  put  into  combining 
the  above  two  methods  in  a  composite  model,  treating  small-scale  roughness  by  the 
Rayleigh-Rice  method  and  large-scale  roughness  by  the  Kirchhoff  method  [3].  A  good 
review  of  theoretical  work  and  experimental  investigations  on  rough  surface  scattering 
can  be  found  in  Ogilvy’s  book  [20]. 

It  is  generally  accepted  that  backscattering  from  sediment  volume  inhomogeneities 
can  contribute  significantly  to  seafloor  backscattering,  especially  at  low  frequencies, 
small  grazing  angles  and  where  the  bottom  is  flat.  Prom  the  1960’s,  a  number  of  mod¬ 
els  have  been  developed  in  order  to  examine  the  underlying  scattering  mechanisms 
and  to  predict  the  strength  of  the  backscattered  signals. 

Stockhausen  [21]  derived  a  volumetric  backscattering  strength  expression  assum¬ 
ing  that  the  water-sediment  interface  is  flat  and  refracting  (with  the  consequent  crit¬ 
ical  angle  effect),  and  with  the  homogeneous  sediment  containing  a  uniform  set  of 
solid  spherical  particles  which  act  to  scatter  the  acoustic  energy.  Treating  the  small 
spheres  as  uncorrelated  point  scatterers,  he  employed  Morse’s  expression  [22]  which 
is  valid  for  scattering  from  spheres  much  smaller  than  a  wavelength.  In  his  model, 
Stockhausen  represented  all  the  scattering  processes  by  a  single  volume  backscatter¬ 
ing  cross  section  without  further  exploring  any  physical  mechanisms. 

Nolle  et  al.  [23]  at  almost  the  same  time  developed  their  own  model.  They  de- 
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scribed  the  scattering  from  the  sediment  volume  based  on  a  random  distribution  of 
scattering  amplitude  per  unit  volume.  The  scattering  strength  distribution  was  con¬ 
sidered  to  be  centered  at  a  mean  value  with  random  deviations.  For  simplicity,  the 
scattering  autocorrelation  function  was  assumed  to  have  the  form  of  an  exponential 
decay.  The  autocorrelation  length  was  proportional  to  the  particle  size  in  the  sed¬ 
iment  volume,  thus  going  one  step  further  compared  to  Stockhausen’s  uncorrelated 
point  scatterer  model.  Nolle  et  al.  also  conducted  an  experiment  in  the  laboratory 
studying  acoustic  scattering  from  a  smooth  sand  surface  and  compared  their  model 
with  the  collected  data.  It  showed  that  the  model  had  some  difficulties  in  explaining 
the  scattering  from  sub-critical  angles. 

Crowther  [24]  included  both  interface  roughness  and  volume  inhomogeneity  effects 
in  his  ocean  bottom  backscattering  model.  Kuo’s  formula  [25]  for  backscattering  from 
the  rough  interface  between  two  homogeneous  fluids  was  used  here.  Extending  Nolle 
et  a/.’s  model,  he  assumed  the  correlation  function  for  impedance  fluctuations  to  be 
elliptic  exponential,  which  was  able  to  account  for  the  anisotropy  of  the  scatterers. 
However,  this  model  also  had  problems  with  scattering  at  sub-critical  angles  when 
comparison  with  Nolle  et  a/.’s  laboratory  backscattering  data  was  made.  Without 
knowledge  of  the  detailed  shape  of  the  correlation  function,  it  would  be  inappropriate 
to  predict  such  characteristics  as  the  frequency  dependence  of  the  scattering  employ¬ 
ing  this  model. 

Morse  and  Ingard  [17],  studied  volume  scattering  due  to  compressibility  and  den¬ 
sity  fluctuations  in  a  free-space  scenario.  Their  approach  is  still  one  of  the  best 
available  methods  in  modeling  volume  scattering. 

There  were  some  singular  features  of  the  backscattering  coefficient  obtained  in 
experimental  measurements  at  sea  [26]  [27].  These  included  the  frequency  indepen¬ 
dence  (or  weak  frequency  dependence)  over  the  frequency  range  1-100  KHz  and  an 
angular  dependence  proportional  to  sin9  for  grazing  angles  6  from  5  to  50  degrees. 
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In  order  to  interpret  this,  Ivakin  and  Lysanov  [28]  proposed  their  statistical  mod¬ 
els  on  the  basis  of  the  geoacoustical  model  presented  earlier  by  Lysanov  [29].  The 
scattering  was  thought  to  come  from  “sharply  anisotropic  random  inhomogeneities 
(fluctuations  of  the  refractive  index):  large-scale  in  the  horizontal  plane  and  small- 
scale  in  depth”  in  the  sediment.  They  used  the  Born  approximation  to  derive  an 
expression  for  the  equivalent  surface  scattering  coefficient,  which  was  not  dependent 
on  the  frequency  if  the  absorption  coefficient  was  proportional  to  frequency  and  the 
power  spectrum  of  the  inhomogeneities  was  inverse  proportional  to  the  wavenumber 
to  the  third  order.  Also,  they  extended  the  single-scale  (horizontal  correlation  coeffi¬ 
cient)  model  to  two-scale  and  multi-scale  models  so  as  to  account  for  the  variability 
of  the  horizontal  correlation  coefficient  in  different  regions.  The  interface  effect  was 
ignored  by  assuming  that  the  changes  in  sound  speed  and  density  from  the  water  to 
the  water-saturated  sediment  were  very  small.  The  model  prediction  agreed  very  well 
with  the  data  after  the  model/data  fit  of  a  free  parameter:  the  ratio  of  the  vertical 
and  horizontal  correlation  lengths. 

In  their  next  paper  [30],  Ivakin  and  Lysanov  revised  their  model  to  account  for 
the  interface  roughness  effect.  The  Kirchhoff  approximation  was  applied  to  the  rough 
interface  backscattering  case.  The  model  studied  in  detail  the  influence  of  interface 
roughness  on  the  volume  backscattering  cross  section:  for  slow  or  nonreflecting  bot¬ 
toms,  the  volume  scattering  did  not  depend  on  whether  the  bottom  interface  was 
smooth  or  not;  for  fast  bottoms,  the  roughness  effect  was  significant  in  the  sub- 
critical  grazing  angle  region.  Physically,  the  reason  for  the  latter  result  was  because 
the  rough  interface  would  enhance  the  penetration  of  sound  into  the  bottom  medium 
at  small  grazing  angles  and  therefore  intensify  the  scattering  of  sound  by  volume  in¬ 
homogeneities.  They  emphasized  the  nonadditiveness  of  the  scattering  effects  due  to 
volume  inhomogeneities  and  water-bottom  interface  roughness,  which  was  different 
from  several  other  researchers  who  assumed  that  these  two  effects  were  uncorrelated. 
The  authors  did  not  include  the  possible  lateral  wave  effects. 
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Ivakin  [15]  extended  his  model  for  scattering  from  volume  inhomogeneities  again 
to  deal  with  a  stratified  bottom.  The  correlation  function  for  the  random  spatial 
fluctuations  of  the  acoustical  parameters  was  considered  to  be  “statistically  homo¬ 
geneous  in  the  horizontal  plane  and  quasi-homogeneous  with  respect  to  the  vertical 
coordinate”.  The  model  allowed  considerable  changes  in  the  sound  speed  and  the 
density  with  depth.  He  studied  the  linear  increase  in  sound  speed  and  density  case 
and  the  linear  decrease  in  sound  speed  case.  The  results  showed  good  model/data 
agreement. 

For  the  high-frequency  (10-100  kHz)  bottom  backscattering  model  proposed  by 
Jackson  et  al.  [3]  in  1986,  the  composite  roughness  approximation  was  applied  for  the 
scattering  due  to  interface  roughness.  To  include  the  volume  scattering  contribution, 
they  employed  Stockhausen’s  formula  [21]  and  accounted  for  the  volume  scattering  by 
an  equivalent  surface  backscattering  coefficient.  The  assumption  inherent  in  Jackson 
et  al.’s  approach  was  the  neglect  of  any  correlation  between  the  scattering  due  to  in¬ 
terface  roughness  and  sediment  inhomogeneities,  which  was  different  from  Ivakin  and 
Lysanov’s  model  [30].  Multiple  scattering  was  also  ignored.  The  comparison  of  the 
model  with  two  sets  of  data  suggested  that  “in  soft  sediment,  sediment  volume  scat¬ 
tering  is  likely  to  be  more  important  than  roughness  scattering,  except  near  normal 
incidence  and  for  grazing  angles  smaller  than  the  critical  angle.  For  sand  bottoms, 
roughness  scattering  is  relatively  more  important” .  However,  the  volume  scattering 
parameter  was  still  obtained  from  model/data  fits  without  any  relationship  to  the 
sediment  properties. 

Mourad  and  Jackson  [31]  generalized  Jackson  et  a/.’s  model  [3]  by  including  the 
sediment  sound  absorption  in  the  interface  boundary  condition.  They  constructed 
some  empirical  relationships  for  estimating  surficial  values  of  the  input  geoacous¬ 
tic  parameters  to  the  model  using  bulk  measurements  of  logarithmic  grain  size 
proposed  by  Hamilton  [32].  These  parameters  otherwise  could  only  be  obtained  by 
direct  measurement,  which  was  not  available  for  many  datasets.  Nevertheless,  as  they 
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pointed  out,  the  gradients  of  sound  speed  and  density,  which  might  play  an  impor¬ 
tant  role  in  affecting  the  acoustical  response  of  the  sediment,  were  not  included  in 
this  model. 

Later,  Mourad  and  Jackson  [1]  went  one  step  further  by  including  the  gradient 
of  sound  speed  in  their  low-frequency  (100-1000  Hz)  model.  The  volume  scattering 
model  was  very  similar  to  Ivakin’s  model  [15]  and  they  assumed  “a  distribution  of 
uncorrelated  omnidirectional  point  scatterers  in  the  sediment  causes  the  backscatter- 
ing  of  sound”.  The  model  related  the  oscillation  of  the  backscattering  strength  to 
the  acoustic  field  within  the  sediment  and  suggested  that  the  total  backscattering 
from  the  ocean  bottom  was  controlled  by  the  processes  affecting  sediment  volume 
scattering,  for  example  strong  layering.  They  also  discussed  the  possible  errors  in 
the  measurement  of  bottom  scattering  strength  near  the  normal  incidence  direction 
using  omnidirectional  sources  and  receivers.  For  low-frequency  backscattering,  the 
questions  for  this  model  were  how  good  the  local  plane-wave  assumption  inherent 
in  the  model  was  and  the  validity  of  attributing  the  volume  scattering  effect  to  an 
equivalent  surface  scattering  process. 

Hines  [14]  developed  a  backscattering  model  in  an  approach  similar  to  Ivakin’s 
[15].  He  followed  Chernov’s  work  [33]  and  applied  the  Born  approximation  and  the 
far-field  assumption.  Backscattering  of  the  acoustic  energy  in  his  model  was  due  to 
the  sound  speed  and  density  fluctuation  in  the  sediment,  which  he  related  to  the 
fluctuation  of  porosity.  The  lateral  wave  effect  was  for  the  first  time  included  in  this 
model,  attempting  to  interpret  the  phenomena  observed  in  experiments  at  small  graz¬ 
ing  angles.  An  exponential  decay  correlation  function  was  assumed  and  the  model  fit 
several  published  data  very  well.  Yet  a  priori  knowledge  of  the  correlation  function 
and  variance  of  the  porosity  were  needed  for  the  model  prediction.  The  approach  of 
decomposing  the  incident  spherical  wave  into  a  refracting  plane  wave  and  an  evanes¬ 
cent  wave  was  debatable.  Later,  he  extended  the  model  to  deal  with  bistatic  scattering 
[34]. 
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Tang  in  his  thesis  [35]  and  Tang  and  Frisk  in  their  several  papers  [36]  [37]  [38] 
discussed  in  detail  the  scattering  from  a  random  layer  or  half-space  where  the  sound 
speed  is  assumed  to  be  a  constant  plus  a  small  random  component.  A  self-consistent 
volume  scattering  model  using  a  perturbation  method  was  developed  to  model  scat¬ 
tering  from  sound  speed  fluctuations.  A  coherent  field  reflection  coefficient  which 
included  scattering  loss  could  be  calculated.  An  interesting  point  was  that  the  spatial 
correlation  length  of  the  scattered  field  could  be  used  to  infer  the  correlation  length 
of  the  scatterers.  This  provided  a  way  of  inverting  for  the  bottom  parameters  critical 
to  bottom  scattering  by  measuring  the  scattered  field  using  multiple  receivers.  Also 
taken  into  consideration  was  the  anisotropy  of  the  scatterers.  Besides  examining  the 
combination  of  interface  roughness  and  volume  inhomogeneity  effects,  an  attempt  was 
made  to  solve  the  near-field  problem  in  low-frequency  scattering  when  the  far-field 
assumption  was  not  appropriate  anymore. 

Lyons  et  al.  [5]  extended  Jackson  et  aVs  model  [3]  for  scattering  from  the  seafloor. 
In  addition  to  the  composite  roughness  model  for  interface  scattering  and  Stock¬ 
hausen’s  expression  for  volume  scattering  [21],  they  included  the  volume  scattering 
from  a  random  inhomogeneous  continuum  and  scattering  from  subbottom  interfaces. 
Their  approach  to  calculating  the  scattering  from  the  random  inhomogeneous  con¬ 
tinuum  was  similar  to  Hines’  work  [14].  The  compressibility  and  density  variations 
were  modeled.  The  correlation  function  could  have  different  correlation  lengths  in 
the  horizontal  and  in  the  vertical  to  allow  for  anisotropy  in  the  volume  scattering 
model.  They  employed  the  Born  approximation  to  obtain  the  volume  backscattering 
cross  section,  which  was  a  free  parameter  in  Jackson  et  aVs  model,  before  fitting 
it  to  Stockhausen’s  formula  for  an  effective  surface  backscattering  coefficient.  They 
estimated  all  of  the  input  parameters  except  the  horizontal  correlation  length  from 
core  data.  The  comparison  with  data  was  very  good.  Because  this  model  was  more 
complicated  and  had  more  input  parameters,  its  performance  depended  very  much 
on  the  estimation  of  those  parameters. 
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Recently,  Yamamoto  [39]  used  the  same  approach  as  Chernov’s  [33]  to  model 
density  and  sound  speed  fluctuations  in  the  sediment.  The  above  fluctuations  were 
described  by  a  power  law  distribution  and  the  parameters  were  estimated  from  cross¬ 
well  tomography.  The  model  compared  with  data  favorably.  However,  the  ray  tracing 
type  of  propagation  model  would  not  be  able  to  deal  with  scattering  at  low  grazing 
angles  for  a  fast  bottom. 

The  modeling  of  volume  scattering  in  a  shallow  water  waveguide  can  be  found  in 
the  work  by  Ellis  [40],  Tang  [8]  and  Tracey  [7]. 

For  elastic  media,  in  addition  to  compressional  waves,  shear  waves  and  evanescent 
waves  (e.g.,  Scholte  waves  at  a  fluid-elastic  interface)  may  play  important  roles  in 
bottom  scattering.  Their  interactions  generate  some  unique  scattering  features  that 
cannot  be  observed  in  a  fluid-fluid  environment  [41]  [42]  [43]  [44].  For  discrete  elastic 
targets  in  a  layered  waveguide,  the  resultant  scattering  is  treated  by  the  transition 
matrix  formalism  [45]  [46]  or  the  boundary  element  method  (BEM)  [47].  However, 
the  elasticity  effect  will  not  be  considered  in  this  thesis  work. 

On  the  experimental  side,  early  investigations  of  bottom  backscattering  used  an 
omnidirectional  source  and  receiver.  Mackenzie  [48]  presented  the  first  deep-sea  mea¬ 
surement  results  in  1961.  The  source  and  receiver  were  close  to  the  sea  surface.  He 
found  that  the  bottom  backscattering  strength  obeyed  Lambert’s  law  of  diffuse  reflec¬ 
tion  for  grazing  angles  from  30  to  90  degrees,  with  the  scattering  constant  (Mackenzie 
coefficient)  having  a  value  of  -27dB.  This  value  didn’t  change  for  the  frequencies  530 
Hz  and  1030  Hz.  In  1968,  Merklinger  [2]  reported  his  experiment  using  a  hydrophone 
and  explosive  charges  suspended  near  the  ocean  bottom.  The  data  analysis  indicated 
that  the  reverberation  due  to  the  subbottom  layer  structure  contributed  significantly 
to  the  total  bottom  reverberation.  This  highlighted  the  importance  of  subbottom 
inhomogeneities  in  the  prediction  of  the  bottom  backscattered  field. 
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Thanks  to  the  evolution  of  technology,  later  experiments  have  used  directional 
sources  and/or  receivers  to  obtain  scattering  data  from  the  ocean  bottom.  Boehme 
et  al.  [49]  conducted  an  experiment  to  study  backscattering  at  low  grazing  angles 
in  the  frequency  range  30-95  kHz,  and  found  that  the  Lambert’s  law  fit  for  grazing 
angles  ranging  from  2  to  10  degrees.  Preston  and  Akal  [50]  employed  a  towed  hor¬ 
izontal  array  and  a  suspended  vertical  array  to  measure  ocean  basin  reverberation. 
They  found  that  strong  backscattering  regions  were  related  to  basin  topographic  fea¬ 
tures.  Hines  and  Barry  [51]  carried  out  an  experiment  in  the  Sohm  Abyssal  Plain 
with  the  source  and  receiver  suspended  close  to  the  smooth  seabed.  Their  acoustic 
array  contained  an  omnidirectional  hydrophone  and  a  ring  projector  with  vertical 
directivity.  Different  features  were  observed  for  backscattering  from  different  bottom 
types.  Jackson  and  Briggs  [4]  used  a  towed  platform  equipped  with  planar  arrays 
to  study  high  frequency  bottom  backscattering  in  three  different  regions.  Sediment 
volume  scattering  was  found  dominant  in  two  of  the  three  sites. 


1.5  Contributions 

This  thesis  has  advanced  the  state-of-the-art  in  several  ways.  These  include: 

•  Evaluation  of  sound  propagation  models  used  in  bottom  volume  scattering  stud¬ 
ies.  The  validity  and  limitations  of  some  widely  used  models  have  been  recog¬ 
nized  and  the  equivalent  surface  scattering  strength  is  found  insufficient  to  char¬ 
acterize  the  mid-  and  low-frequency  volume  scattering  process  when  multi-path 
effects  are  severe. 

•  Implementation  of  a  numerical  monostatic  backscattering  model  with  the  ability 
to  incorporate  3-D  volume  inhomogeneities  and  to  compare  with  experimental 
data  effectively.  Density  fluctuations  have  been  confirmed  to  be  the  main  con¬ 
tributor  to  backscattering,  while  the  anisotropy  of  the  statistical  distribution  of 
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the  inhomogeneities  in  both  horizontal  and  vertical  directions  is  found  to  have 
a  profound  influence  on  the  angular  dependence  of  the  backscattering  strength. 
In  addition,  the  synthesized  time  series  can  help  to  design  an  experiment  and 
provide  predictions  of  reverberation  fields. 

•  Formulation  of  the  spectral  method  for  the  generation  of  azimuthally  integrated 
3-D  random  fields.  This  approach  provides  the  opportunity  to  understand  the 
full  3-D  phenomena  with  a  manageable  computational  load. 

•  Comparisons  with  the  ARSRP  backscattering  data.  The  good  fit  between  the 
model  and  data  has  proved  the  versatility  of  the  model.  The  power  law  type  of 
power  spectrum  is  found  to  be  the  best  description  of  sound  speed  and  density 
fluctuations  in  the  sediment  that  gives  rise  to  excellent  model/data  comparisons. 

1.6  Overview  of  the  Thesis 

Chapter  2  begins  by  reviewing  the  definition  of  the  scattering  cross  section  and  the 
propagation  models  by  Hines,  Ivakin  and  Mourad  and  Jackson.  Comparisons  with 
the  exact  solution  obtained  by  OASES  are  then  carried  out  for  three  different  type 
of  fluid  bottom  models.  The  pros  and  cons  of  each  propagation  model  are  discussed. 
Meanwhile,  the  validity  of  the  equivalent  surface  scattering  strength  in  characterizing 
volume  scattering  is  also  evaluated. 

In  Chapter  3,  the  formulations  for  scattering  due  to  sound  speed  and  density  fluc¬ 
tuations  are  derived.  In  addition,  the  power  spectra  and  correlation  functions  that 
are  going  to  be  investigated  later  are  presented.  Analytic  solutions  for  scattering 
in  a  free-space  scenario  are  used  to  study  the  eflfects  on  backscattering  imposed  by 
different  distributions  and  different  correlation  lengths  of  sound  speed  variations. 

Chapter  4  presents  a  scheme  to  generate  3-D  azimuthally-summed  random  fields 
by  the  spectral  method.  Numerical  examples  are  also  given  for  verification  of  this 
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approach. 


Chapter  5  describes  the  procedures  for  generating  backscattered  time  series.  Then 
the  focus  shifts  to  the  verification  of  the  numerical  model.  The  effect  of  density  fluc¬ 
tuations  is  studied  at  the  end. 

In  Chapter  6,  we  present  a  detailed  description  of  the  ARSRP  backscattering  ex¬ 
periment.  A  brief  introduction  of  the  data  processing  methodology  is  then  given. 
The  estimation  of  backscattering  strength  from  data  and  the  selection  of  model  pa¬ 
rameters  are  then  discussed.  Data  from  the  ARSRP  experiment  are  compared  to 
numerical  simulations. 

Chapter  7  summarizes  the  physical  insights  gained  from  the  modeling  and  sug¬ 
gests  directions  for  future  work. 
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Chapter  2 


Evaluation  of  Sound  Propagation 
Models  Used  in  Bottom  Volume 
Scattering  Studies 


The  proper  evaluation  of  sound  propagation  between  sources/receivers  and  scatter- 
ers  is  important  in  characterizing  bottom  volume  scattering.  In  this  chapter,  several 
sound  propagation  models  used  in  bottom  volume  scattering  studies  are  evaluated  and 
their  results  compared  to  the  exact  solution  obtained  through  a  numerical  wavenum¬ 
ber  integration  technique.  It  is  found  that  Hines’  approach  [14]  works  well  for  the 
two  isovelocity  half-space  case  except  when  the  grazing  angle  is  close  to  the  critical 
angle.  The  far-field  approximation,  given  by  Ivakin  [15]  and  Mourad  and  Jackson  [1], 
has  a  performance  depending  upon  the  sound  speed  structure  in  the  sediment.  For 
an  isovelocity  slow  bottom,  it  agrees  well  with  the  exact  solution.  However,  discrep¬ 
ancies  arise  for  an  isovelocity  fast  bottom  or  a  bottom  with  a  complex  sound  speed 
structure.  In  addition,  the  appropriateness  of  using  the  equivalent  surface  scattering 
strength  as  a  function  of  grazing  angle  in  volume  scattering  characterizations  is  stud¬ 
ied.  In  conclusion,  precautions  need  to  be  taken  in  modeling  both  the  propagation 
effects  and  the  scattering  mechanisms  associated  with  the  bottom  volume  scattering 
process. 
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This  chapter  is  structured  as  follows:  In  the  first  section,  we  discuss  the  classi¬ 
cal  scattering  problem  and  the  definition  of  the  conventional  scattering  cross  section. 
The  concept  of  an  equivalent  surface  scattering  strength  which  incorporates  bottom 
volume  scattering  is  also  reviewed;  for  the  second  section,  we  describe  the  propaga¬ 
tion  models,  and  in  the  next  section,  we  compare  the  results  of  these  models  for  some 
typical  cases  of  interest. 


2.1  Statement  of  the  Problem 

Acoustic  wave  scattering  from  the  ocean  bottom  generally  includes  rough  water /- 
bottom  interface  scattering  and  subbottom  volume  scattering.  While  rough  surface 
scattering  has  been  the  focus  for  many  years,  recent  evidence  shows  that  volume  scat¬ 
tering,  due  to  inhomogeneities  and/or  scattering  layers  within  the  sediment,  could 
contribute  significantly  to  bottom  reverberation. 

The  conventional  quantities  used  to  characterize  the  above  two  processes  are  the 
surface  scattering  coefficient  and  the  volume  scattering  coefficient  [52],  which  histori¬ 
cally  are  defined  in  the  context  of  a  plane  wave  being  scattered  by  scatterers  confined 
in  an  otherwise  homogeneous  medium.  In  ocean  bottom  scattering,  however,  the 
incident  field  is  normally  not  a  plane  wave.  The  multipath  and  refractive  effects 
resulting  from  bottom  sound  speed  structure  complicate  the  scattering  process  and 
make  it  extremely  difficult  to  pinpoint  the  scattering  element  in  the  classical  sense. 
This  problem  will  not  diminish  when  the  bottom  volume  scattering  process  is  treated 
by  an  equivalent  surface  scattering  process  with  the  parameter  of  equivalent  surface 
scattering  coefficient  being  used.  Therefore,  the  credibility  of  the  scattering  strength 
estimated  from  experimental  data  by  this  means  needs  to  be  examined. 

Bottom  volume  scattering  modeling  usually  include  two  components:  the  scatter¬ 
ing  part,  which  encompasses  the  appropriate  scattering  mechanisms  and  the  propa- 
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gation  part,  which  takes  into  account  propagation  from  the  source  to  the  scatterer 
and  from  the  scatterer  to  the  receiver.  We  will  concentrate  on  propagation  mod¬ 
els  in  this  study,  which  is  equivalent  to  evaluating  the  Green’s  function  between  the 
source/receiver  and  the  scatterer.  The  necessity  to  include  poro-elastic  effects  in  the 
scattering  model  remains  an  open  question  [18],  and  we  will  regard  the  bottom  as  an 
acoustic  fluid.  The  free-space  Green’s  function,  together  with  the  transmission  coeffi¬ 
cient  at  the  water/bottom  interface  has  been  used  in  some  work  [24,  53].  Despite  its 
simplicity,  it  cannot  incorporate  the  interface  wave  contribution  at  subcritical  grazing 
angle  for  a  fast  bottom.  Since  we  are  interested  only  in  the  methods  that  are  good 
for  the  entire  angular  regime,  this  method  will  not  be  discussed  in  this  paper.  Ivakin 
[15]  and  Mourad  and  Jackson  [1]  considered  the  incident  wave  on  the  interface  to  be 
a  plane  wave,  with  the  incident  direction  varying  over  the  insonified  region  to  match 
the  true  incident  direction.  For  convenience,  we  call  it  the  far-field  approximation 
in  this  paper.  This  approach  is  actually  a  simplified  version  of  the  stationary  phase 
method  with  only  the  most  dominant  stationary  point  contribution  being  included. 
Meanwhile,  Hines  [14]  followed  the  stationary  phase  approach  of  Westwood  [54]  to 
find  the  transmitted  field  due  to  a  point  source  in  the  two  isovelocity  half-space  case. 

The  goal  of  this  chapter  is  to  demonstrate  the  validity  and  limitations  of  the 
above-mentioned  propagation  models  based  on  realistic  geoacoustic  parameters.  Re¬ 
sults  from  these  models  are  compared  with  the  exact  solution  obtained  through  a 
numerical  wavenumber  integration  technique.  Their  influence  on  the  estimation  of 
the  equivalent  surface  scattering  strength  is  also  studied  in  this  chapter. 


2.2  The  Scattering  Cross  Section 

The  conventional  scattering  cross  section  is  defined  when  describing  the  scattering 
phenomenon  associated  with  a  single  scatterer  (e.g.,  a  particle),  as  shown  in  Fig.  2-1. 
For  simplicity,  the  incident  plane  wave  is  assumed  to  be  propagating  in  the  positive 
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Figure  2-1;  Geometry  of  the  classical  scattering  problem. 


z  direction  with  unit  amplitude  and  suppressed  time  dependence 

(2.1) 

The  time  dependence  is  chosen  in  the  same  convention  as  that  in  Ivakin’s  and  Hines’ 
work.  However,  bear  in  mind  that  I  will  change  it  to  starting  from  Chapter  3  in 
order  to  be  consistent  with  the  convention  used  in  statistic  theory  and  random  field 
generation.  The  scattered  wave  behaves  as  a  spherical  wave  in  the  far  field  and  is 
measured  at  infinity: 

gtfcr 

=  forr-^oo,  (2.2) 

where  r,  6  and  cf)  are  spherical  coordinates  with  the  origin  at  the  scatterer.  The 
scattering  cross  section  is  then  defined  as  [55] 

<^{^,(1))  =  (2.3) 

It  is  a  measure  of  the  scattered  power  in  the  {d,  4>)  direction  per  unit  solid  angle,  per 
unit  incident  intensity.  Because  the  incident  field  is  a  plane  wave  and  the  observa- 
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As  stated  before,  in  bottom  scattering  studies,  scattering  by  the  rough  water/- 
bottom  interface  and  volume  inhomogeneities  within  the  bottom  are  two  key  pro¬ 
cesses.  Any  type  of  scattering  modeling  will  inevitably  encounter  the  problem  of 
describing  the  scattering  ability  of  a  surface  area  and/or  a  volume.  Figure  2-2  shows 
the  concept  underlying  the  definitions  of  a  surface  scattering  coefficient  Og  and  a  vol¬ 
ume  scattering  coefficient  <t„,  which  is  similar  to  o.  Plane  wave  incidence  is  again 
assumed  here,  with  the  source  and  receiver  in  the  far  field.  The  incident  intensity  /q  is 
measured  at  unit  distance  away  from  the  scatterer.  The  scattered  wave  intensity  Ig  is 
obtained  at  the  receiver  and  then  multiplied  by  in  oder  to  compensate  for  spherical 
spreading  of  the  scattered  wave.  Similar  to  the  definition  of  cr,  the  scattering  ability 
is  characterized  by  the  ratio  of  the  scattered  wave  intensity  to  the  incident  intensity 
at  unit  distance  from  the  scatterer  per  unit  area  or  volume.  It  can  be  expressed  in 
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the  following  formula  [52]: 


(Js  =  Is{0A)r‘^/[h{0(iAQ)dS],  (2.4) 

and 

Gy  =  Is(e,  (f>)r^/[Io{9o,  (j)o)dV],  (2.5) 

where  ($0,  ^o)  is  the  incident  direction,  (6,  <j))  is  the  scattering  direction,  and  dS  and 
dV  are  the  insonified  area  and  volume,  respectively. 

Several  things  need  to  be  noticed  about  the  introduction  of  Gs  and  Gy.  First, 
instead  of  a  single  scatterer,  there  exists  a  distribution  of  scatterers  in  the  area  S  or 
volume  V .  Therefore,  Gs  and  Gy  represent  a  measurement  of  the  average  scattering 
ability  of  the  insonified  surface  or  volume.  Second,  in  the  definition  of  a,  a  plane  wave 
incident  on  the  scatterer  is  assumed  in  order  to  exclude  the  effect  of  source  position. 
Similarly  for  Gg  and  Gy,  the  incident  wavefront  is  considered  planar  or  locally  planar, 
with  the  source  in  the  far  field.  However,  due  to  the  possible  multipath  and  refractive 
effects,  this  is  generally  not  true  for  bottom  scattering  problem.  Third,  for  the  same 
reasons,  the  scattered  wave  is  not  simply  a  spherical  wave.  As  a  result,  the  spreading 
loss  at  a  distance  r  cannot  be  compensated  for  by  the  factor  r^,  and  in  general,  the 
phase  is  incorrectly  calculated  as  well. 

Despite  the  increasing  interest  in  bottom  volume  scattering,  it  remains  a  diffi¬ 
cult  issue  to  characterize  the  scattering  process.  Unlike  the  rough  surface  scatterers 
which  are  concentrated  at  the  water/bottom  interface,  the  scatterers  associated  with 
volumetric  scattering  may  be  distributed  over  a  large  region.  Scatterers  may  be  lo¬ 
cated  immediately  beneath  the  water/bottom  interface  or  a  certain  distance  below 
the  seafloor.  As  illustrated  in  Fig.  2-5,  if  we  have  a  fast  bottom  and  the  scatterer  is 
not  at  the  interface  but  with  jd  smaller  than  critical  grazing  angle,  the  scattered  wave 
can  reach  the  receiver  through  two  paths,  therefore  at  two  different  angles.  The  path 
lengths  are  different  so  that  the  arrival  times  are  not  the  same.  As  will  be  discussed 
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Figure  2-3;  Schematic  illustration  of  the  equivalent  surface  scattering  coefficient  def¬ 
inition. 

later,  the  above  multi-path  effect  will  pose  nontrivial  problems  for  scattering  strength 
estimations. 

When  the  bottom  attenuation  is  sufficiently  high  and  the  depth  of  acoustic  pen¬ 
etration  is  much  smaller  than  the  distance  between  the  source  and  the  sediment, 
scattering  from  volume  inhomogeneities  within  the  bottom  can  be  described  by  an 
equivalent  surface  scattering  process  [1].  As  shown  in  Fig.  2-3,  scattering  from  the 
shaded  region,  which  is  a  slice  of  the  scattering  volume  in  an  axisymmetric  coordi¬ 
nate  system,  is  attributed  to  scattering  from  the  surface  area  above  it.  All  of  the 
scattering  are  considered  all  from  the  grazing  angle  9  because  the  difference  would  be 
minimal  when  the  observation  point  (i.e.,  receiver)  is  very  far  away.  The  equivalent 
surface  scattering  coefficient  agg,  instead  of  the  volume  scattering  coefficient  a„,  is 
then  defined  to  quantify  the  scattering  process.  Mourad  and  Jackson  [1]  present  a 
thorough  description  of  this  approach,  with  the  result 

aes  =  J  ayl'^il^dz,  (2.6) 
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where 


'J'l  =  (2.7) 

Here  2:  is  the  depth  coordinate,  and  are  the  incident  field  as  it  would  be  mea¬ 
sured  at  the  volume  element  dV  and  the  area  element  dA,  respectively,  a^s  can  be 
viewed  as  a  measure  of  the  scattering  contribution  from  the  whole  shaded  region. 
is  the  normalized  field  as  if  the  scatterers  were  at  the  surface. 

The  concept  of  an  equivalent  surface  scattering  coefficient  has  been  used  exten¬ 
sively  in  scattering  data  analysis.  In  the  high-frequency,  high-attenuation  case,  the 
incident  wave  may  penetrate  to  relatively  shallow  depths.  However,  at  low  frequen¬ 
cies  and  attenuations,  the  ensonification  can  encompass  a  substantial  portion  of  the 
sediment  column.  The  sediment  layering  and  the  sound  speed  gradient  cause  the 
scatterers  to  be  ensonified  by  an  incident  field  consisting  of  refracted,  multiple  re¬ 
flected  paths.  Even  if  we  neglect  the  multiple  scattering  effect,  which  is  to  ignore  the 
rescattering  of  the  scattered  wave,  it  is  still  very  difficult  to  relate  the  scattered  re¬ 
turns  recorded  in  an  experimental  time  series  to  the  corresponding  scatterers.  As  can 
be  seen  in  Fig.  2-3,  the  thickness  of  the  shaded  region  would  be  very  large  for  deep 
penetration.  In  the  above  calculation  of  the  equivalent  surface  scattering  strength, 
all  the  contributions  from  each  shaded  “vertical  bar”  would  be  ascribed  to  only  one 
scattering  angle,  which  can  be  justified  for  shallow  penetration  in  a  high-frequency 
situation  but  is  invalid  in  the  low-frequency  case.  In  fact,  the  volume  element  could 
affect  the  scattering  at  different  angles  due  to  multipath  effects.  As  a  result,  the 
angular  dependence  of  the  predicted  scattering  strength  might  be  error-prone.  In 
a  shallow  water  environment,  the  limited  water  depth  would  further  deteriorate  the 
problem,  since  the  distance  between  the  source/receiver  and  the  water/bottom  inter¬ 
face  is  confined.  All  of  the  above  would  directly  affect  the  effectiveness  of  using  the 
scattering  strength  as  a  criteria  in  a  model/data  comparison. 
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2.3  Bottom  Propagation  Models 


Propagation,  i.e.,  the  Green’s  function,  plays  an  important  role  in  the  modeling  of 
the  scattering  process.  In  this  section,  three  methods  will  be  described  for  calculating 
the  propagation  from  a  point  source  in  the  water  column  to  a  scatterer  in  the  sedi¬ 
ment.  The  propagation  from  the  scatterer  to  a  receiver  in  the  water  can  be  obtained 
accordingly  by  the  principle  of  reciprocity. 

The  wave  equation  in  an  inhomogeneous  medium  [56]  can  be  written  as 

"W  *>)  -  =  /(>•.  <).  (2-8) 

where  p(r,  t)  represents  the  pressure,  p(r)  represents  the  density,  c(r)  represents  the 
sound  speed,  and  /(r,  t)  represents  the  source  term  as  a  function  of  space  r  and  time 
t.  If  constant  density  and  harmonic  time  dependence  [exp(— iwi)]  are  assumed,  Eq.2.8 
becomes  the  Helmholtz  equation 

[V^  +  A:o(r)(l  -I-  e(r))]  P(r,  w)  =  F(r, a;),  (2.9) 

where  A;o(r)  is  the  medium  wavenumber  for  the  background  sound  speed  Co(r)  and 
frequency  oj,  with 

Mr)  =  (2.10) 

and 


p(r,t)  =  P(r,a;)e--*. 

Parameter  e(r)  is  the  perturbation  to  the  wavenumber  due  to  the  inhomogeneous 
sound  speed.  Using  a  perturbation  method,  the  total  field  P(r)  (for  simplicity,  a;  will 
be  implicit  from  now  on)  can  be  expressed  as  the  summation  of  the  incident  field 
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Po(r)  and  the  scattered  field  Ps(r): 


P(r)  =  Po(r)  +  P,(r),  (2.11) 

where  the  incident  field  satisfies 

[V^  +  kl(r)]  Po(r)  =  P(r).  (2.12) 

We  assume  a  point  source  at  r  =  and  the  Green’s  function  can  be  obtained  from 

[V^  +  A:^(r)]  ^(r,  Fs)  =  5(r  -  Fs),  (2.13) 

and  the  boundary  conditions.  The  scattered  field  P^  satisfies 

[V^  +  A:g(r)]  P,(r)  =  -A:^(F)e(F)P(F).  (2.14) 

The  integral  equation  for  the  solution  of  the  scattered  field  is 

P,(r)  =  -A:^(f)  j  G(f,  F')e(F')P(r')dr'.  (2.15) 


From  Eq.2.15,  it  can  be  seen  that  the  scattered  field  and  the  Green’s  function  are 
closely  related.  The  Green’s  function,  which  characterizes  the  propagation  from  the 
source  to  the  scatterer  and  from  the  scatterer  to  the  receiver,  plays  an  important  role 
in  modeling  of  the  bottom  scattering  process.  In  this  section,  three  methods  will  be 
described  for  calculating  the  propagation  from  a  point  source  in  the  water  column  to 
a  scatterer  in  the  sediment.  The  propagation  from  the  scatterer  to  a  receiver  in  the 
water  can  be  obtained  accordingly  by  the  principle  of  reciprocity. 

There  are  several  ways  to  solve  Eq.2.13.  For  simplicity,  we  will  consider  a  range- 
independent,  horizontally  stratified  environment,  which  If  we  choose  a  cylindrical 
coordinate  system,  with  the  vertical  2:-axis  passing  through  the  source  and  the  r-axis 
being  parallel  to  the  interfaces,  it  becomes  an  axisymmetric  propagation  problem.  At 


46 


Figure  2-4:  The  experimental  scenario. 

this  point,  we  will  assume  unit  source  strength.  The  Green’s  function  can  then  be 
calculated  by  means  of  wavenumber  integration.  For  the  scenario  shown  in  Fig.  2-4, 
we  have  the  Green’s  function  in  the  sediment  G(r,  z)  as  follows, 

fOO 

G{r,z)  =  /  g{kr,z)krJo{krr)dkj. 

Jo 

=  ^  /  9{kr,z)krH^^\krr)dkr  (2.16) 

Z  J—oo 

where  kr  is  the  horizontal  wavenumber,  Jq  is  the  zeroth  order  Bessel  function,  is 
the  zeroth  order  Hankel  function  of  the  first  kind  and  g{kr,z)  satisfies  the  equation 

^  4-  {k‘^{z)  -  k^)^g{kr,z)  =  0  (2.17) 

and  the  boundary  conditions  of  continuity  of  normal  particle  velocity  and  pressure 
across  each  interface.  When  krV  ':$>  1,  we  have 

1  too  f— 

=  9{kr,z)y/kre^^'‘^^-*'>dkr  (2.18) 

y/zirr  J-oo 

Eq.2.18  can  be  evaluated  numerically.  The  procedure  of  implementation  can  be 
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found  in  Schmidt’s  SAFARI  manual  [16]  and  Ref. [57].  In  fact,  the  numerical  solution 
of  the  wave  equation  in  Section  III  is  obtained  primarily  by  using  SAFARI.  One  of 
the  advantages  of  SAFARI  is  that  it  generally  provides  stable  results  for  complicated 
sediment  structures. 


A  point  source  is  located  at  (0, 2:5).  If  an  isovelocity  water  column  is  assumed  so 
that  the  water  wavenumber  ko  =  oj/cq,  the  depth-dependent  Green’s  function  in  the 
water  column  will  be 

where  k^  is  the  vertical  wavenumber  in  the  water.  By  matching  the 

boundary  conditions  at  the  water /sediment  interface  2;  =  0,  it  is  clear  that  the 
depth-dependent  Green’s  function  in  the  sediment  can  be  expressed  in  terms  of 
_gifczol^d/47ri/j^^.  In  other  words,  we  have 


gifciQ  \zi\ 

g{kr,  z)  =  -gi{kr,  ^  , 


(2.20) 


where  giikrik)  can  be  viewed  as  the  solution  to  Eq.2.17  with  a  plane  wave  incident 
on  the  water/sediment  interface.  Using  the  coordinates  shown  in  Fig.  2-4,  Eq.2.18 
now  becomes 


Gir,  z)  =  -=  /  ; - dkr-  (2.21) 

V27rr  J-00  ^T^kzo  ^ 

The  stationary  phase  method  is  often  used  to  evaluate  an  integral  of  this  type. 

Following  Ivakin’s  approach  [15],  one  obtains 


G{t,z)  = 


(2.22) 


where  R  is  the  distance  between  the  source  and  the  point  at  the  interface  right  above 
the  scatterer  in  the  sediment,  as  shown  in  Fig.  2-5.  For  simplicity,  a  two  half-space 
scenario  is  being  considered.  In  the  above  approach,  the  propagation  can  be  inter¬ 
preted  as  a  spherical  wave  traveling  from  the  source  to  the  interface  point  B  above 
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Figure  2-5:  Schematic  of  scattering  geometry  for  a  two  isovelocity  half-space  case. 

the  scatterer  and  then  a  plane  wave  traveling  from  point  B  to  the  scatterer.  For  a 
slow  bottom,  there  is  only  one  eigenray  path.  Therefore,  Path  OBS  is  used  to  approx¬ 
imate  Path  OA  ’S.  Point  B  is  the  approximate  stationary  point  and  A  ’  is  the  exact 
one;  For  a  fast  bottom,  there  are  two  cases.  When  the  incident  grazing  angle  a  is 
larger  than  the  critical  angle,  the  situation  is  similar  to  that  for  a  slow  bottom.  Only 
one  eigenray  path  exists  and  A  is  the  exact  stationary  point.  However,  when  a  is 
smaller  than  the  critical  angle,  both  Path  OAS  and  OBS  will  exist.  Paths  OAS  and 
OBS  correspond  to  the  refracted  wave  and  the  the  evanescent  wave,  respectively.  If 
cx  is  close  to  the  critical  angle,  there  might  be  more  than  two  eigenrays  according  to 
Westwood  [54].  On  the  one  hand,  point  B is  not  the  same  as  point  A  or  A’,  which  will 
bring  in  a  certain  degree  of  error  for  a  larger  than  the  critical  angle.  On  the  other 
hand,  for  a  larger  than  the  critical  grazing  angle,  the  stationary  point  A  contribution, 
i.e.,  the  refracted  wave,  will  be  neglected.  When  the  bottom  attenuation  is  high  and 
the  scatterer  is  close  to  the  water /sediment  interface,  this  approach  is  justified,  since 
the  ray  path  of  the  refracted  wave  will  be  much  longer  than  that  of  the  lateral  wave. 
The  refracted  wave  would  be  greatly  attenuated  because  of  that.  However,  we  will 
be  able  to  see  some  differences  if  the  above  conditions  cannot  be  satisfied.  Mourad 
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and  Jackson’s  approach  [1]  is  somewhat  similar  to  Ivakin’s.  Instead  of  regarding  the 
propagation  as  consisting  of  a  spherical  wave  traveling  from  the  source  to  the  inter¬ 
face  point,  the  incident  wave  is  considered  to  be  a  plane  or  locally  plane  wave  with 
incident  direction  allowed  to  vary  over  the  insonified  region  in  order  to  match  the 
true  incident  direction.  For  convenience,  we  will  call  the  above  approach  the  far-field 
approximation  in  this  work. 

The  formula  that  Hines  derived  in  his  paper  [14],  also  a  far-field  approach,  works 
for  a  two  isovelocity  half-space  environment  only.  He  applied  the  stationary  phase 
method  as  well.  In  order  to  distinguish  it  from  the  above-mentioned  far-field  approx¬ 
imation,  we  will  name  it  Hines’  method  in  the  later  comparison.  In  a  two  isovelocity 
half-space  situation,  gi(kr,z)  in  Eq.2.21  can  be  further  expressed  as 

^i(A:r,2)  =  (2.23) 

where  T^s  is  the  transmission  coefficient  from  the  water  to  the  sediment,  is  the 
vertical  wavenumber  in  the  sediment,  and  2  is  the  distance  from  the  interface  to  the 
scatterer.  If  Eq.2.21  is  rewritten  in  terms  of  the  angles  as  plotted  in  Fig.  2-5  and 
includes  the  result  in  Eq.2.18,  we  will  have 

■ir  I  kn  /■— too+7r/2  , -  - - 

G{r,  z)  =  e*  4  y  (2.24) 

where  Oi  is  the  incident  angle  and  =  ki/ko-  h  here  represents  the  medium 
wavenumber  in  the  sediment.  Eq.  2.24  is  a  little  different  from  that  in  Hines’  paper 
since  we  choose  Oi  to  be  the  incident  angle  from  the  water  to  the  sediment  instead  of 
the  opposite.  The  phase  of  the  integrand  is  now 

/ (Oi)  =  iko{rcosei  -  z^sinOi  -f  z^ri^^  -  cos^Oi).  (2.25) 

n^s  here  can  be  complex  to  account  for  the  attenuation  in  the  sediment.  According 
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to  the  stationary  phase  method,  the  stationary  point  can  be  found  from  the  following 


dm 

dOi 


=  —iko  i^sinOi  +  ZgCosOi  — 


zsin{29i) 


(2.26) 


There  will  be  two  solutions  to  Eq.2.26  for  a  fast  bottom.  The  solution  corre¬ 
sponding  to  the  evanescent  wave  can  be  obtained  by  perturbing  the  incident  angle  13, 
where 

P  =  tan~^(—Zs/r).  (2.27) 


The  real  solution  can  be  found  by  performing  a  first-order  Taylor  expansion, 


de^p- 


f'iP) 

fipy 


(2.28) 


where  the  prime  (')  represents  the  derivative  with  respect  to  the  incident  angle.  We 
have  to  use  a  complex  root  finder  to  obtain  the  second  solution  6r,  which  corresponds 
to  the  refracted  wave.  The  perturbation  method  suggested  by  Hines  has  some  difiScul- 
ties  when  the  actual  stationary  point  is  close  to  a  singular  point,  the  critical  angle  in 
this  case.  It  might  also  be  a  problem  for  9e  and  is  possibly  a  cause  of  the  overestimate 
for  the  field  in  the  region  close  to  the  critical  angle  as  can  be  seen  later. 


For  the  two  half-space,  fast  bottom  situation,  the  incident  spherical  wave  is  there¬ 
fore  decomposed  into  two  plane  waves;  the  refracted  wave  and  the  evanescent  wave. 
The  evanescent  wave  exists  only  in  the  subcritical  grazing  angle  region.  The  station¬ 
ary  phase  method  gives  the  results 

(2.29) 

where 


A{9)  =  Ty,s  (9)  y/^e^ko{rcos9-ZsSine+zy/nl, -cos^d)  ^ 


(2.30) 
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and 


B(e)  = 


N 


r  (  rcos9  —  z.sinO  —  z 


cos(26) 


sin^{2d) 


cos26>  AyJ (n2  ^  -  cos'^ey 


).  (2.31) 


For  the  refracted  wave 

Gr  =  QiOr), 


(2.32) 


and  for  the  evanescent  wave 

Ge^Qiee). 


(2.33) 


Prom  the  results  in  the  next  section,  this  method  can  be  seen  to  be  reasonably 
good  for  the  two  isovelocity  half-space  case,  except  in  the  region  close  to  the  critical 
angle,  where  there  are  more  than  two  stationary  points,  as  pointed  out  by  Westwood 
[54].  However,  it  is  not  trivial  to  find  all  the  stationary  points  for  a  bottom  with  mul¬ 
tiple  layers  or  with  a  sound  speed  gradient  because  the  phase  term  of  the  integrand 
in  Eq.2.19  becomes  much  more  complicated. 


2.4  Comparison  of  Propagation  Models 

The  far-field  approximation  and  Hines’  method  are  widely  implemented  in  bottom 
scattering  modeling,  especially  in  high-frequency  problems.  Yet  their  validity  and 
limitations  need  to  be  recognized  when  dealing  with  general  experimental  scenarios. 
In  this  section,  based  on  real  geoacoustic  parameters,  the  above  two  propagation  mod¬ 
els  will  be  evaluated  through  comparisons  with  exact  solutions  obtained  by  numerical 
wavenumber  integrations.  An  effort  will  also  be  made  to  identify  their  possible  influ¬ 
ence  on  the  quantification  of  bottom  scattering. 

We  will  concentrate  on  three  types  of  fluid  models  [cf.  Fig.  2-6]:  two  isovelocity 
half-spaces;  two  half-spaces  with  an  upward-refracting  sound  speed  profile  in  the  sed¬ 
iment;  a  layer  with  an  upward-refracting  profile  between  two  isovelocity  half-spaces. 
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Figure  2-6:  The  sound  speed  profile  for  three  types  of  bottom:  (a)  two  isovelocity 
half  spaces;  (b)  two  half  spaces  with  a  sound  speed  gradient  in  the  bottom;  (c)  three 
layers  with  a  sound  speed  gradient  in  the  middle  layer. 
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Most  of  the  geoacoustic  parameters  used  in  the  comparisons  are  measurement  results 
or  extrapolations  of  measured  data  which  have  appeared  in  the  literature.  In  this  sec¬ 
tion,  /  will  stand  for  frequency,  a  and  p,  will  stand  for  the  sound  speed  and  density 
in  the  overlying  water  {i  =  0),  just  beneath  the  water /sediment  interface  {i  =  1),  and 
in  the  lower  half  space  (i  =  2)  in  Fig.  2-6(c).  The  density  is  chosen  to  be  constant 
in  each  layer.  The  loss  factor  5i,{i  =  0,1,2)  is  the  ratio  of  the  imaginary  and  real 
parts  of  the  complex  wavenumber  at  the  top  of  each  sediment  layer.  The  quantity  g 
represents  the  sound  speed  gradient,  and  H  stands  for  the  middle  layer  thickness  in 
the  three-layer  case.  The  wavelength  in  the  water  is  given  by  A,  and  2:^  is  the  source 
height  above  the  water/bottom  interface  and  z  is  the  distance  between  the  scatterer 
in  the  bottom  and  the  water/sediment  interface.  The  quantity  Zs  and  z  are  scaled 
relative  to  wavelength  A.  The  source  height  z^  is  chosen  to  be  lOOA  in  order  to  avoid 
the  near-field  effect,  which  is  true  for  the  whole  section.  The  transmission  loss  in  the 
comparisons  is  defined  as  the  ratio  in  decibels  between  the  acoustic  pressure  at  the 
scatterer  for  a  simple  point  source  and  the  pressure  produced  at  a  distance  of  Im 
from  the  same  source. 


2.4.1  Two  Isovelocity  Half-Spaces 

The  sound  speed  in  the  sediment  depends  on  the  sediment  type.  The  sediment /water 
sound  speed  ratio  could  be  larger  or  smaller  than  1,  which  determines  the  existence 
of  a  critical  angle.  According  to  Hamilton  [58],  for  soft  sediments  such  as  silty  clay, 
the  sound  speed  just  beneath  the  water/sediment  interface  is  usually  less  than  that  of 
the  overlying  water.  The  geoacoustic  parameters  used  in  Fig.  2-7,  which  is  consistent 
with  Continental  Rise  data  and  taken  from  Fig.  3  in  Mourad  and  Jackson’s  paper 
[1],  is  an  example.  From  the  field  at  scatterers  situated  at  z  =  0  and  z  =  -2A,  it 
can  be  seen  that  results  of  the  far-field  approximation  and  Hines’  method  are  almost 
identical  to  that  of  the  numerical  wavenumber  integration.  There  is  no  critical  angle 
and  only  one  stationary  point,  i.e.,  scatterers  are  insonified  by  the  refracted  wave 
only.  Both  approaches  are  good  approximations  in  this  situation.  One  would  expect 
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Figure  2-7:  Comparison  of  transmission  loss  in  the  environment  with  no  critical 
angle.  Parameters  are:  /  =  500  Hz,  Cq  =  1530m/s,  Ci  =  1510m/s,  pi  =  1.72g/cm^, 
6  =  0.00164,  Zg  =  lOOA.  The  transmitted  field  at  (a)  z  =  0;  (h)  z  =  — 2A. 

the  calculation  of  the  equivalent  surface  scattering  strength  to  be  appropriate  and  ac¬ 
curate  without  including  the  multipath  effect.  Since  Zg  and  z  are  both  scaled  relative 
to  A,  the  results  obtained  here  and  for  the  two  isovelocity  half-space  case  should  be 
frequency-independent  if  the  bottom  attenuation  is  linear  with  frequency,  which  is  an 
open  question  [59,  60,  32]. 

For  sandy  bottoms  and  turbidite  bottoms  in  certain  areas,  the  sediment /water 
sound  speed  ratio  is  larger  than  1.  Input  parameters  in  Fig.  2-8  correspond  to  the 
Hudson  Canyon  experiment  reported  by  Rogers  et  al.  [61],  where  the  bottom  is  de¬ 
scribed  as  silty  sand  and  silty  clay.  We  choose  the  frequency  to  be  l^kHz  and  assign 
the  bottom  attenuation  to  be  O.ldB/m/kHz  ((5i  =  0.0027),  which  is  at  the  high  end 
for  this  type  of  bottom.  The  transmitted  field  is  again  shown  at  different  depths  in 
Fig.  2-8.  It  is  found  that  the  agreement  between  Hines’  method  and  the  numerical 
wavenumber  integration  technique  is  fairly  good  except  for  some  overestimate  in  the 
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Figure  2-8:  Comparison  of  transmission  loss  with  the  critical  angle  effect  for  a  low 
attenuation  bottom.  Parameters  are:  /  =  500  Hz,  cq  =  1500m/s,  ci  =  1560m/s, 
Pi  =  1.8^/cm^,  S  =  0.0027,  Zg  =  lOOA.  The  transmitted  field  at  (a)  z  =  0;  (b) 
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Figure  2-9:  The  detailed  comparison.  Parameters  are  the  same  as  that  in  Fig.  2-8(c). 

region  just  below  the  critical  grazing  angle.  The  critical  angle  {9c  =  cos“^(co/ei))  is 
16  degrees  and  is  indicated  in  the  figure.  The  far-field  approximation,  on  the  other 
hand,  tends  to  underestimate  the  field  at  subcritical  grazing  angles.  Differences  will 
increase  with  the  depth  of  the  scatterer  [cf.  Fig.  2-8].  In  Fig.  2-9,  we  add  the 
evanescent  wave  and  refracted  wave  contributions  to  Fig.  2-8  (c).  For  the  far-field 
approximation,  it  is  clear  that  the  neglect  of  the  refracted  wave  in  the  subcritical 
angle  region  is  the  cause  of  the  underestimate.  The  small  error  before  the  critical 
angle  can  be  attributed  to  the  difference  between  the  true  stationary  point  and  that 
used  in  the  approximation.  As  for  Hines’  method,  the  match  is  excellent  everywhere 
except  in  a  small  region  after  the  critical  angle.  This  is  the  transition  area  from  a  field 
composed  of  one  stationary  point  to  a  field  composed  of  two  stationary  points,  where 
the  approximation  would  normally  be  invalid,  according  to  Westwood  [54].  Another 
possible  cause  might  be  the  failure  of  the  first-order  perturbation  method  to  find  the 
stationary  point  close  to  the  singularity  associated  with  the  critical  angle. 
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Figure  2-10:  The  critical  angle  effect  for  a  high  attenuation  bottom.  Parameters 
are:  /  =  40000  Hz,  cq  =  1500  m/s,  Ci  =  1689m/s,  pi  =  1.97g/cm^,  S  =  0.0166, 
Zs  =  lOOA.  The  transmitted  field  at  (a)  ^  =  0;  (b)  ^  =  -A;  (c)  2  =  -2A. 


For  a  sandy  bottom,  the  attenuation  will  be  much  higher  than  that  for  a  soft  bot¬ 
tom.  The  geoacoustic  parameters  in  Fig.  2-10  are  from  the  CBBL  Panama  City  site 
[18].  Frequency  of  AQkHz  is  chosen,  the  same  as  that  in  the  experiment.  It  is  a  sandy 
bottom  and  the  attenuation  was  found  to  be  {i.h2AdB/m/kHz,  which  is  equivalent  to 
5i  —  0.0166  .  The  critical  angle  is  now  27.4  degrees.  The  figure  shows  results  similar 
to  those  in  Fig.  2-8  for  the  field  at  three  different  depths.  Nevertheless,  the  differ¬ 
ence  between  the  far-field  approximation  and  the  numerical  wavenumber  integration 
is  smaller  than  that  in  the  lower  attenuation  bottom.  The  interference  between  the 
refracted  wave  and  the  evanescent  wave  is  hardly  noticeable  at  subcritical  angles. 
The  explanation  for  this  is  that  the  refracted  wave  has  a  longer  ray  path  than  that 
of  the  evanescent  wave  and  is  therefore  greatly  attenuated  at  long  ranges  in  the  high- 
attenuation  bottom.  In  other  words,  the  evanescent  wave  contribution  will  prevail  at 
subcritical  angles.  So  we  see  from  the  above  that  the  single  stationary  point  treat¬ 
ment  in  the  far-field  approximation  works  better  in  a  high-attenuation  bottom  than 
in  a  low-attenuation  one. 

The  effect  of  different  propagation  models  on  the  estimation  of  scattering  strength 
estimation  is  of  considerable  interest  in  bottom  scattering  modeling.  In  Figs.  2-11 
and  2-12,  we  have  calculated  the  equivalent  surface  scattering  strength  [cf.  Eq.2.6] 
associated  with  volume  backscattering  using  the  same  parameters  as  those  in  Figs. 
2-8  and  2-10.  Point  scatterers  are  assumed  here  for  simplicity,  and  the  volume  scat¬ 
tering  cross  section  cr„  is  chosen  to  be  0.00003  in  Fig.  2-11  and  0.000906  in  Fig.  2-12. 
Both  in  the  ARSRP  low-frequency  sediment  scattering  experiment  and  the  CBBL 
high-frequency  scattering  experiment  at  the  Eckernfoerde  Bay  site,  strong  scattering 
layers  are  found  to  be  a  certain  distance  away  from  the  water/sediment  interface 
[9,  6].  So  we  will  choose  a  scattering  layer  with  a  thickness  of  50A  here  starting  from 
different  depths  below  the  interface.  Figures  2-11,  2-12(b)(d)(f)  show  the  difference 
in  backscattering  strength  between  the  exact  solution  and  the  far-field  approximation 
result.  Since  results  from  Hines’  method  fit  exact  solutions  very  well,  but  in  a  very 
narrow  region,  we  consider  that  it  would  yield  similar  outputs  and  will  not  discuss  it 
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Figure  2-11:  The  equivalent  surface  backscattering  strength  for  a  low  attenuation 
bottom  assuming  point  scatterers.  Parameters  are  the  same  as  that  in  Fig.  2-8  and 
=  0.00003.  The  backscattering  strength  and  the  difference  as  the  scattering  layer 
begins  at  (a)(b)  z  =  0;  (c)(d)  z  =  -A;  (e)(f)  z  =  -2A. 
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Figure  2-12:  The  equivalent  surface  backscattering  strength  for  a  high  attenuation 
bottom  assuming  point  scatterers.  Parameters  are  the  same  as  that  in  Fig.  2-10  and 
ay  =  0.000906.  The  backscattering  strength  and  the  difference  as  the  scattering  layer 
begins  at  (a)(b)  z  =  0;  (c)(d)  z  =  -A;  (e)(f)  z  =  -2 A. 


here. 


From  Fig.  2-11,  it  is  evident  that  the  underestimate  of  the  field  at  subcriti- 
cal  angles  for  the  far-field  approximation  leads  to  differences  in  equivalent  surface 
backscattering  strength.  However,  the  biggest  difference  is  only  about  IdB  when  the 
scattering  layer  starts  right  at  the  interface.  On  the  other  hand,  when  the  scatter¬ 
ing  layer  begins  at  greater  depths,  the  difference  is  larger.  For  the  scattering  layer 
starting  at  two  wavelengths  beneath  the  interface,  up  to  IhdB  difference  appears. 
Also,  we  notice  that  the  actual  level  of  the  backscattering  strength  in  this  case  is 
only  several  dB  lower  than  that  for  the  scattering  layer  starting  from  the  interface. 
Therefore  one  can  hardly  distinguish  them  from  the  scattering  measurement  without 
a  priori  knowledge,  which  highlights  the  necessary  precautions  in  scattering  modeling 
using  the  far-field  approximation  in  this  environment.  When  the  bottom  attenuation 
is  low,  the  penetration  is  deep.  The  multipath  effect  is  inevitable.  The  scattered 
returns  will  reach  the  receiver  via  different  paths,  therefore  different  scattering  an¬ 
gles,  which  makes  the  angular  dependence  of  the  scattering  strength  less  meaningful. 
All  of  these  results  cast  some  doubt  on  the  validity  of  using  the  equivalent  surface 
scattering  strength  to  characterize  the  volume  scattering  process. 

Fig.  2-12  shows  similar  results  for  a  sandy  bottom.  The  difference  in  backscatter¬ 
ing  strength  is  smaller  when  compared  with  that  in  Fig.  2-11.  Numerical  errors  can 
even  be  seen  in  Fig.  2-12(e)  and  (f).  The  scattering  level  for  a  scattering  layer  deep 
in  the  sediment  is  much  lower  than  that  at  the  interface.  Since  the  sand  is  usually 
a  strong  scatterer  and  the  bottom  penetration  is  much  shallower  due  to  high  atten¬ 
uation,  the  problems  that  are  addressed  regarding  Fig.  2-11  are  much  less  serious. 
Both  the  far-field  approximation  and  the  equivalent  surface  scattering  strength  are 
appropriate. 
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2.4.2  Two  Half-Spaces  with  Sound  Speed  Gradient  in  the 


Lower  Half 

Sound  speed  gradients  exist  in  both  the  water  column  and  sediment  in  real  ocean 
environments  [62]  and  have  significant  effects  on  sound  propagation.  The  gradient  in 
the  water  column  is  generally  small  compared  with  that  in  the  sediment  and  will  be 
neglected.  We  will  therefore  focus  on  the  effect  of  the  sound  speed  gradient  in  the 
sediment.  We  will  not  include  Hines’  method  in  the  comparisons  from  now  on  due  to 
the  complexity  of  extending  it  to  complicated  sound  speed  structures. 

In  order  to  obtain  an  analytic  solution,  we  choose  the  l/(?  linear  sound  speed 
profile, 

Ci(z)  =  ci(l  —  (2.34) 

for  a  real  constant  a,  with  units  of  inverse  length.  The  sound  speed  gradient  at  the 
interface  is 

dci{z)  cio 

^  dz  z=o  2 

Figure  2-13  shows  the  frequency  dependence  of  the  comparison  between  the  far- 
field  approximation  and  the  numerical  wavenumber  integration  for  the  case  where 
there  is  a  drop  in  sound  speed  at  the  water/sediment  interface.  The  sound  speed 
gradient  g  has  a  typical  value  1.0s“^.  The  mismatch  of  the  oscillation  pattern  is 
obvious  at  lOOHz,  but  as  the  frequency  increases,  the  match  is  better.  At  500IIz, 
they  are  almost  indistinguishable.  The  oscillation  here  comes  from  the  interference 
between  the  direct  incident  wave  and  the  upward-refracted  wave.  The  sound  speed 
gradient  will  determine  the  path  length  of  the  upward  refracting  ray.  At  high  fre¬ 
quencies,  the  refracted  wave  will  be  attenuated  so  much  that  virtually  no  interference 
occurs.  Therefore  the  higher  the  frequency,  the  better  the  far-field  approximation. 
The  elucidation  of  the  mismatch  in  oscillation  patterns  rests  on  the  examination  of 
the  refracted  wave.  In  the  far-field  approximation,  the  field  is  dominated  by  the  in¬ 
terference  between  the  plane  wave  incident  at  a  particular  angle  and  the  refracted 
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Figure  2-13:  The  frequency  dependence  of  the  comparison  results.  Parameters 
Co  =  1530m/s,  Cl  =  1510m/s,  pi  =  l.72glcm^,  5  =  0.00164,  g  =  ls-\  =  1 
2  =  0  and  (a)  /  =  100  Hz-,  (b)  /  =  300  Hz-,  (c)  /  =  500  Hz. 


wave  originating  from  the  same  angle  of  incidence.  While  in  the  numerical  wavenum¬ 
ber  integration  approach,  the  interference  involves  refracted  waves  originating  from 
many  incident  angles.  It  all  boils  down  to  the  plane  wave  assumption  in  the  far-field 
approximation,  which  is  not  included  in  the  wavenumber  integration  approach. 

The  parameters  in  Fig.  2-14  correspond  to  the  case  where  there  is  an  increase  in 
sound  speed  at  the  water/sediment  interface.  The  gradient  g  is  again  chosen  to  be 
1.0s“^.  Fig.  14  shows  the  field  for  three  different  frequencies  at  two  different  depths. 
As  we  anticipated,  the  field  now  includes  the  contribution  from  the  evanescent  wave, 
the  interface  refracted  wave  and  the  turning  refracted  wave,  as  shown  in  Fig.  2-15. 
Together,  they  will  determine  the  interference  pattern.  The  critical  angle  effect  is  still 
clear  [cf.  Sec.  2.4.1].  At  higher  frequencies,  the  results  approach  those  for  the  two 
isovelocity  half-space  case  because  the  turning  refracted  wave  is  strongly  attenuated. 
All  in  all,  the  far-field  approximation  will  have  some  problems  at  low  frequency  for 
this  kind  of  low-attenuation  bottom  with  a  sound  speed  gradient.  As  for  the  high 
attenuation  sandy  bottom,  we  would  expect  similar  phenomena  as  those  at  high  fre¬ 
quency  in  Fig.  2-14,  as  can  be  seen  in  Fig.  2-16.  The  parameters  in  Fig.  2-16  are  the 
same  as  in  for  Fig.  2-10  except  that  the  frequency  is  lOOi/z.  Also,  we  have  studied 
the  effect  of  the  sound  speed  gradient.  The  gradient  will  decide  the  path  length  for 
the  turning  refracted  wave,  therefore  the  amplitude,  in  a  bottom  with  attenuation. 
We  would  expect  higher  amplitude  for  the  turning  refracted  wave  in  a  bottom  with  a 
larger  gradient,  because  the  path  length  would  be  shorter.  As  a  result,  the  oscillation 
is  more  significant  for  a  larger  gradient  as  shown  in  Fig.  2-17. 

Again  we  will  study  the  influence  of  the  propagation  model  on  the  quantification 
of  the  scattering  process.  We  still  assume  point  scatterers  here  and  calculate  only 
the  backscattering  strength.  The  scattering  layer  will  be  considered  starting  from  the 
water /sediment  interface  only.  For  a  high-attenuation  sandy  bottom  like  that  used 
to  generate  Fig.  2-16,  we  choose  to  be  0.000906.  As  shown  in  Fig.  2-18,  the 
equivalent  surface  backscattering  strength  calculated  with  the  far-field  approxima- 
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Figure  2-14.  The  frequency  and  depth  dependence  of  the  comparison  results  for  a 
low  attenuation  bottom  with  a  critical  angle.  Parameters  are  cq  =  1500m/s,  ci  = 
1560m/s,  Pi  =  1.8p/cm3,  6  =  0.0027,  g  =  ls-\  =  lOOA  and  (a)  /  =  100  Hz,z  =  0 
;  (b)  /  =  100  Hz,  z  =  -2A;  (c)  f  =  bOO  Hz,  z  =  0  ■,  (d)  /  =  500  Hz,  z  =  -2A;  (e) 
/  =  5000  Hz,z  =  Q-,{i)  f  =  5000  Hz,  z  =  -2A. 
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Figure  2-15:  The  propagation  paths  for  two  half-space  fast  bottom  with  an  upward- 
refracting  sound  speed  profile  in  the  bottom. 


Figure  2-16:  The  depth  dependence  for  a  high-attenuation  bottom.  Parameters  are: 
/  =  100  Hz,  Co  =  1500  m/s,  Ci  =  1689m/s,  pi  =  1.97p/cm^,  5  =  0.0166,  g  =  ls“^ 
2:5  =  lOOA.  The  transmitted  field  at  (a)  2  =  0;  (b)  z  =  — 2A. 
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Figure  2-17:  The  sound  speed  gradient  dependence.  Parameters  are  the  same  as  that 
in  Fig.  2-13(a)  except  the  gradient  (a)  g  =  (b)  g  =  4s~^;  (c)  g  =  10s“^ 


Figure  2-18:  The  equivalent  surface  backscattering  strength  for  a  high-attenuation 
bottom  assuming  point  scatterers.  Parameters  are  the  same  as  that  in  Fig.  2-16 
with  the  starting  depth  of  the  scattering  layer  to  be  at  z  =  0.  (a)  the  backscattering 
strength;  (b)  differences. 
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Figure  2-19:  The  equivalent  surface  backscattering  strength  for  a  low-attenuation 
bottom  assuming  point  scatterers.  Parameters  are:  /  =  2300  Hz,  cq  =  1550  m/s, 
Cl  =  1576m/s,  Pi  =  1.83g/cm\  6  =  0.00193,  g  =  lOs-^,  z,  =  lOOA.  The  scattering 
layer  starts  at  z  =  0.  (a)  the  backscattering  strength;  (b)  differences. 

tion  is  very  close  to  the  exact  solution.  For  a  low-attenuation,  soft  sediment  bottom, 
we  let  the  geoacoustic  input  be  the  same  as  that  for  Fig.  15  (b)  in  Mourad  and 
Jackson’s  paper  [1].  Fig.  2-19  shows  the  backscattering  strength  predicted  by  the 
far-field  approximation  and  the  numerical  wavenumber  integration  technique.  Notice 
that  the  oscillations  are  absent  for  the  numerical  solution  but  still  present  for  the  far- 
field  approximation.  An  interesting  phenomenon  appearing  in  the  results  obtained 
by  numerical  wavenumber  integration  is  that  the  backscattering  strength  has  a  drop 
of  about  3dB  at  16  degrees,  which  is  not  shown  in  the  results  when  the  far-field  ap¬ 
proximation  is  applied.  In  Mourad  and  Jackson’s  figure  15(b),  the  experimental  data 
also  shows  an  8  -  lOdB  fall.  Although  the  drop  in  our  prediction  does  not  match  the 
data  exactly,  it  is  one  step  closer  and  suggests  that  the  far-field  approximation  may 
not  be  good  enough  in  situations  like  this. 
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Figure  2-20:  Figure  15(b)  in  Mourad  and  Jackson’s  paper  [1].  The  circles  here  repre¬ 
sent  the  data  points. 
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2.4.3  Three  Layers  with  a  Sound  Speed  Gradient  in  the  Mid¬ 
dle  Layer 

A  more  realistic  model  of  the  sediment  is  to  have  a  basement  beneath  the  wa¬ 
ter/sediment  interface.  We  still  choose  a  1/c^  linear  sound  speed  profile  below  the 
water /sediment  interface  but  constant  sound  velocity  in  the  basement.  Scatterers  are 
assumed  to  be  distributed  in  the  middle  layer.  The  presence  of  an  additional  layer 
interface  enables  the  insonification  of  scatterers  through  some  extra  paths,  i.e.,  the 
single  or  multiple  bounces  between  the  two  interfaces,  which  is  similar  to  some  degree 
to  the  situation  where  there  is  a  sound  speed  gradient  in  the  bottom.  The  turning 
point  there  can  be  considered  to  be  the  virtual  lower  boundary. 

One  would  gather  that  the  middle  layer  thickness  H  is  likely  to  play  an  important 
role  in  the  interference  between  multipath  rays.  Fig.  2-21  shows  the  field  at  the 
water/sediment  interface  for  different  middle  layer  thicknesses.  The  parameters  for 
the  middle  layer  are  the  same  as  those  for  Fig.  2- 13(c)  and  correspond  to  a  slow 
bottom.  For  the  thin  middle  layer,  the  “frequency”  of  the  oscillation  is  lower  than 
that  of  the  thick  layer.  However,  the  agreement  between  the  far-field  approximation 
and  the  exact  solution  is  better.  The  “frequency”  of  the  oscillation  is  related  to  the 
rate  of  the  phase  change  for  the  transmitted  wave  with  the  change  of  the  incident 
angle.  For  a  thicker  middle  layer,  a  small  change  of  the  incident  angle  would  cause 
a  big  change  of  the  phase  of  the  transmitted  then  reflected  wave,  which  results  in 
high  “frequency”  oscillations.  The  interference  pattern  will  depend  on  the  phase  of 
the  reflected  rays  and  therefore  on  the  layer  thickness.  The  geoacoustic  input  for  Fig. 
2-22  for  the  upper  two  layers  is  the  same  as  that  for  Fig.  2-12(b),  corresponding  to 
a  sandy  bottom.  Although  the  argument  for  the  “frequency”  of  the  oscillation  still 
holds,  the  level  of  oscillation  is  much  smaller  because  of  the  high  attenuation.  In 
conclusion,  the  far-field  approximation  works  better  for  a  high-attenuation  than  for 
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Figure  2-21:  The  middle  layer  thickness  dependence  for  a  three  layer  environment 
with  a  low  attenuation  middle  layer.  Parameters  are:  /  =  500  Hz,  Cq  =  1530m/s, 
Cl  =  1510m/s,  C2  =  1800m/s,  pi  =  1.72p/cm^,  p2  =  =  0.00164, 

82  =  0.019,  g  =  ls-\  Zs  =  lOOA,  2  =  0  and  (a)  H  =  2A;  (b)  H  =  10 A. 
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Figure  2-22:  The  middle  layer  thickness  dependence  for  a  three-layer  environment 
with  a  high  attenuation  middle  layer.  Parameters  are:  /  =  100  Hz,  cq  =  1500m/s, 
Cl  =  1689m/s,  C2  =  1800m/s,  pi  =  l-OTp/cm^,  pg  =  2.2glcm^,  5i  =  0.0166,  82  = 
0.019,  g  =  4s-i,  =  lOOA,  z  =  0  and  (a)  H  =  2 A;  (b)  H  =  lOA. 
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a  low-attenuation  bottom  in  three  layer  case. 


2.5  Summary 

In  this  chapter,  we  have  compared  several  propagation  models  that  are  widely  used 
in  scattering  studies  with  the  exact  solution  obtained  by  a  numerical  wavenumber 
integration  technique. 

Hines’  method,  which  decomposes  the  spherical  wave  into  two  plane  waves:  the 
refracted  wave  and  the  evanescent  wave,  matches  the  exact  solution  very  well  in  a  two 
isovelocity  half-space  scenario.  The  overestimate  of  the  field  in  a  small  region  close 
to  the  critical  angle  only  affects  the  results  at  less  than  two  degrees  of  grazing  angle 
in  our  comparison,  which  is  tolerable.  Yet  the  difficulty  in  obtaining  the  stationary 
points  easily  and  accurately  hampers  its  generalization  to  more  complex  sound  speed 
structure. 

The  far-field  approximation,  on  the  other  hand,  can  be  easily  implemented  in 
complicated  bottom  models,  partly  because  it  includes  only  one  stationary  point  and 
has  simplified  its  stationary  point  searching  procedure  by  sacrificing  some  accuracy. 
As  a  result,  the  agreement  with  the  exact  solution  is  somewhat  case-dependent.  In  a 
two  isovelocity  half-space  model,  the  approach  works  almost  perfectly  if  there  doesn’t 
exist  the  critical  angle  at  the  water/sediment  interface,  which  is  likely  to  be  the  case 
for  a  great  amount  of  soft  sediment  types  such  as  mud  and  silty  clay.  However,  in 
some  environments,  the  presence  of  the  critical  angle  for  such  low-attenuation,  soft 
sediment  will  find  some  underestimate  of  the  field  at  subcritical  angles.  The  situation 
is  worse  for  deeper  scatterers.  This  is  because  of  the  omission  of  the  refracted  wave 
past  the  critical  angle  inherent  in  the  far-field  approximation.  The  high  attenuation 
in  sandy  bottoms  would  ease  the  problem  to  some  extent  by  attenuating  the  refracted 
wave  due  to  its  long  path.  As  for  the  influence  on  the  equivalent  surface  scattering 
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strength  calculation,  results  are  encouraging  if  the  scattering  layer  starts  from  the 
water/sediment  interface.  However,  when  the  scattering  layer  starts  some  distance 
away  from  the  interface,  which  is  unlikely  for  a  sandy  bottom  but  possible  for  soft 
sediments,  the  error  introduced  by  applying  the  far-field  approximation  is  significant. 

The  occurrence  of  the  sound  speed  gradient  in  the  sediment  provides  an  additional 
propagation  path,  i.e.,  the  turning  refracted  wave,  which  contributes  to  the  determi¬ 
nation  of  the  interference  pattern.  The  attenuation  together  with  the  frequency  and 
the  gradient  controls  the  amplitude  and  phase  of  the  turning  refracted  wave.  There¬ 
fore,  the  results  of  comparisons  are  frequency-  and  gradient-dependent.  In  general, 
the  far-field  approximation  fares  better  for  a  high-attenuation,  low-gradient  bottom. 
While  for  a  low-attenuation  bottom,  the  difference  is  considerably  large,  as  can  be 
seen  in  the  equivalent  surface  scattering  strength  calculation. 

Sediment  layering  is  common  in  a  real  ocean  environment.  In  the  three  layer  case, 
the  bounce  of  the  acoustic  wave  between  two  interfaces  adds  to  multipath  phenom¬ 
ena.  The  oscillation  of  the  field  is  sensitive  to  the  layer  thickness  no  matter  what 
type  of  sediment  is  in  the  middle  layer.  High  bottom  attenuation  is  again  helpful 
in  alleviating  the  multipath  effect,  i.e.,  the  oscillation.  The  far-field  approximation 
matches  the  numerical  wavenumber  integration  better  for  a  thinner  middle  layer  and 
a  higher  attenuation  situation. 

The  validity  of  using  the  equivalent  surface  scattering  strength  to  characterize  the 
volume  scattering  process  is  debatable.  Because  of  the  evident  multipath  contribu¬ 
tion  to  the  field  in  many  cases  discussed  in  the  comparison,  the  obtained  angular 
dependence  of  the  equivalent  surface  scattering  strength  is  under  question.  However, 
as  shown  in  our  results,  it  is  still  reasonable  to  apply  the  equivalent  surface  scattering 
strength  concept  in  some  high-frequency  and  high-bottom  attenuation  bottom  scat¬ 
tering  studies. 
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Chapter  3 


Formulations  for  Modeling  Bottom 
Volume  Backscattering 


A  sound  wave  can  be  scattered  by  a  region  in  which  the  acoustic  properties  of  the 
medium  differ  from  their  values  in  the  rest  of  the  medium.  Sound  speed  and  den¬ 
sity  fluctuations  are  generally  considered  the  major  contributors  to  volume  scattering. 
Although  the  effect  of  density  fluctuation  is  insignificant  within  the  water  column,  ac¬ 
cording  to  Chernov  [33],  it  is  often  the  dominant  factor  in  the  case  of  bottom  backscat¬ 
tering  from  sediment  volume  inhomogeneities  [14].  Prom  a  practical  point  of  view,  a 
statistical  approach  seems  to  be  the  only  resort  for  people  to  model  such  small  scale 
features.  The  fluctuations  of  sound  speed  and  density  in  the  bottom  are  usually  small 
due  to  the  nature  of  sedimentation  process.  Therefore,  the  method  of  small  pertur¬ 
bation  is  suitable  and  widely  used  in  volume  scattering  modeling  [7, 14,  17,  28,  33,  35]. 

In  this  chapter  we  will  derive  the  formula  for  modeling  volume  scattering  due  to 
sound  speed  and  density  fluctuations.  It  is  similar  to  Morse  and  Ingard’s  [17],  Cher¬ 
nov’s  [33]  and  Hines’  [14]  approaches  and  used  by  Tang  in  the  study  of  shallow  water 
reverberation  [8].  Also,  one  section  will  concentrate  on  various  power  spectral  den¬ 
sity  functions  and  correlation  functions  that  are  going  to  be  used  in  describing  sound 
speed  and  density  fluctuations.  In  another  section,  we  will  present  the  analytical 
solution  for  backscattering  from  different  kinds  of  random  sound  speed  fluctuations 
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in  a  free-space  scenario. 


3.1  Volume  Scattering  Theory 


As  shown  in  the  previous  chapter,  the  wave  equation  in  an  inhomogeneous  medium 
can  be  written  as 


p„V.(^Vp(M)) 


1  aV(r,  t) 
c2(r)  dt^ 


/(r,  t), 


(3.1) 


where  p(r,<)  represents  the  pressure,  p(r)  represents  the  density,  c(r)  represents  the 
sound  speed,  and  / (r,  t)  represents  the  source  term  as  a  function  of  space  r  and  time  t. 
Assuming  a  point  source  and  harmonic  time  dependence  exp{iu}t),  the  inhomogeneous 
Helmholtz  equation  is  [33]; 


[V^  +  k\T)]  P(T)  -  ^  .  VP(r)  =  -  r.), 


(3.2) 


where  the  source  strength  is  S^-  We  assume  that  the  background  density  is  constant 
and  the  background  sound  speed  only  varies  with  depth  With  small  variations  in 
both  density  and  sound  speed,  we  have 


P{r)  =  Po  +  Sp{t) 

c(r)  =  Co(2:) +(5c(r),  (3.3) 

where  6p{t)  <C  po  and  5c(r)  <C  co{z)  and  they  both  have  zero  means. 


If  we  expand  the  Helmholtz  equation  and  keep  up  to  first-order  of  smallness  with 
respect  to  p  and  c,  Eq.3.2  becomes 


(V.  ^  _  mm .  VP(.)  = 

Co[z)  Po 


0- 


(3.4) 
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The  total  pressure  field  can  also  be  expressed  as  the  summation  of  the  mean  field 
Po{t)  and  the  scattered  field  Ps{r): 


P(r)  =  Po(r)  +  P.{r). 


(3.5) 


We  can  obtain  the  equation  for  the  mean  field  Po(r)  by  inserting  this  summation  in 
the  equation  above  and  averaging: 


+  kUz)]Poir)  -  { 


25c{r)kQ{z) 

Co(z) 


n(r)) 


Po 


(3.6) 


Subtracting  this  equation  from  the  full  averaged  Helmholtz  equation,  we  can  get  an 
equation  for  the  scattered  field: 


[V^  +  kl(zmr)  = 


2Sc{T)kQ{z) 

Co{z) 


Po(r)  + 


V(^p(r)) 


Po 


VPo(r). 


(3.7) 


It  is  interesting  to  notice  that  the  two  terms  on  the  RHS  caused  by  sound  speed 
and  density  fluctuation  act  like  virtual  sources.  As  discussed  by  Morse  and  Ingard 
[17],  the  first  term  comes  from  the  fact  that  the  small  region  at  r  doesn’t  have  the 
same  sound  speed  (therefore  compressibility)  as  the  surrounding  medium,  which  re¬ 
sults  in  a  monopole-like  source  distribution.  The  second  term  gives  rise  to  dipole-like 
scattering  phenomenon  because  the  small  region  doesn’t  move  in  response  to  the  force 
VPo  with  the  same  velocity  as  does  the  surrounding  medium  due  to  the  difference 
in  density.  This  dipole  source  would  have  a  null  in  the  forward  scattering  direction 
but  its  major  axis  is  in  the  backscattering  direction,  which  makes  density  fluctuations 
a  significant  contributor  to  backscattering.  As  a  result,  it  is  of  great  importance  to 
include  density  variations  in  volume  backscattering  modeling. 


For  the  weak  volume  scattering  problem,  a  usual  assumption  is  that  the  effect  of 
the  scattered  field  on  the  mean  field  can  be  ignored  (Born  Approximation).  Therefore, 
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Eq.3.6  becomes: 


[V^  +  /c^(2:)]Po(r)  =  SJ{t  -  Ts).  (3.8) 

The  corresponding  Green’s  function  can  be  obtained  from 

[V2  +  A:^(r)]G(r,rs)  =  (5(r-rs),  (3.9) 

and  the  boundary  conditions. 


Also,  the  sound  speed  variation  and  density  variation  are  often  considered  linked, 
which  would  simplify  the  problem.  As  discussed  by  Hines  [14],  for  a  random  sediment 
bottom,  both  fluctuations  can  be  expressed  as  functions  of  the  bottom  porosity  p: 


Sc  =  ^Sp 
dp 

Sp  =  pp, 


where  Sp  is  the  variation  in  bottom  porosity.  If  we  use  the  notation 


(3.10) 


e(r) 

Po 


Sc 

Co 

Hr), 


(3.11) 


where  /5  is  the  ratio  between  density  variation  and  sound  speed  variation,  the  scattered 
field  equation  becomes 


[V"  +  kl{z)]Ps{T)  =  2kUz)e{T)Po{T)  +  /?Ve(r)  •  VPo(r).  (3.12) 

Applying  Green’s  Theorem  to  the  equation  above,  we  have 

Ps(Rr)  =  [2kl{z')e{T')Po(T')  +  /?V€(r')  •  VPo(r')]G(Rr,  r')dr',  (3.13) 
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where  Rp  stands  for  the  receiver  position  and  r'  stands  for  the  scatterer  position. 
After  integration  by  parts,  the  scattered  field  equation  can  be  further  expressed  as 

P,(R,)  =  /  £(r')[(2  +  /3)A;2(2')n(r')G(Rx,r')  -/3(VP„(r')  •  VG(R„r')))*'.(3.14) 

Jv 

The  scattered  intensity  4  can  then  be  obtained  by 

j.  =  <in(R.)i^> 

=  /  /  (£(ri)£'(r2))l(2  +  /J)fcJ(z.)P„(r,)G(R„r,)-/3(VF„(ri).VG{R„r,))] 

Jvi  Jv2 

[(2  +  fi)kl{z,)P;{r:,)G’(Rr,  rj)  -  /^(VPo'CrJ  ■  VG-{R„r2)))*,*j.  (3.15) 

Here  (*)  stands  for  complex  conjugate  operation. 

In  the  above  equation,  (e(ri)e*(r2))  is  the  correlation  function  of  the  sound  speed 
fluctuation.  Since  the  value  of  the  sound  speed  variation  must  be  a  real  number,  (*) 
can  be  omitted  from  the  correlation  function  expression.  In  this  weak  volume  scat¬ 
tering  model  with  the  Born  approximation,  we  can  see  that  the  scattered  intensity 
is  determined  by  the  second-order  statistics  of  the  random  sound  speed  and  density 
variations  if  the  propagation  part  is  known.  So  it  is  worthwhile  to  devote  the  next 
section  to  different  kinds  of  correlation  functions  to  be  used  in  volume  scattering 
modeling. 

If  scattering  time  series  are  of  interest,  we  can  apply  the  technique  of  Fourier 
synthesis.  The  scattered  return  y{t)  can  be  calculated  by 

1  r°° 

y{t)  —  —  Ps{'Rr,(^)S{(jo)exp{iu;t)du,  (3.16) 

zTT  J —oo 

where  uj  is  the  radial  frequency  and  S{uj)  is  the  spectrum  of  the  source  signal.  This 
equation  will  be  used  in  Chapter  5  to  generate  scattering  time  series. 

The  validity  of  the  volume  scattering  formula  is  not  restricted  to  layered  fluid 
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media.  As  far  as  scatterers  are  located  in  the  fluid  layer,  the  bottom  can  be  fluid, 
fluid-elastic,  or  elastic  layered  media.  The  random  scattering  layer  is  not  necessarily 
immediately  beneath  the  water /sediment  interface.  The  attenuation  of  both  com- 
pressional  waves  and  shear  waves  and  the  background  sound  speed  gradient  are  easily 
incorporated  into  the  model  because  they  are  taken  care  of  in  the  calculation  of  the 
Green’s  function  by  OASES.  This  theory  will  be  suitable  for  studying  scattering  by  a 
sediment  bottom,  where  interface  roughness  does  not  necessarily  dominate  the  scat¬ 
tering  process.  This  is  exactly  the  situation  that  we  encountered  in  the  ARSRP 
scattering  experiment. 


3.2  Correlation  Functions  of  Sound  Speed  Varia¬ 
tions 

The  second-order  statistics  of  sound  speed  and  density  variations  are  of  great  impor¬ 
tance  to  volume  scattering  modeling,  as  can  be  seen  in  the  section  above.  However,  it 
is  difficult  to  have  direct  measurement  of  such  sediment  properties,  let  alone  a  good 
model  of  them.  Traditionally,  Gaussian  and  exponential  correlation  functions  have 
often  been  used  because  of  their  mathematical  simplicity  [14,  33,  35].  Recently,  a 
power-law  type  of  correlation  function  has  attracted  more  interest  [53,  63].  Mean¬ 
while,  Yamamoto,  in  his  measurement  of  velocity  variability  using  crosswell  acoustic 
tomography,  found  the  measured  power  spectra  approximate  a  power  law  [64].  Still 
it  is  hardly  conclusive  that  the  power-law  type  power  spectral  density  is  the  best 
description  of  velocity  variations  in  the  sediment.  Therefore,  we  choose  to  include  all 
three  of  them  in  our  model. 

In  the  beginning  of  this  chapter,  we  have  defined  two  random  variables,  the  relative 
sound  speed  variation  e  and  the  density  variation.  Since  they  are  assumed  to  be  fully 
correlated,  only  the  statistics  of  e  will  be  dealt  with  from  now  on.  In  addition  to  the 
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zero  mean  assumption  regarding  e,  we  require  that  the  sound  speed  fluctuation  be 
a  homogeneous  random  variable.  In  other  words,  the  process  is  spatially  wide-sense 
stationary,  by  which  we  mean  that  the  covariance  of  the  random  variable  is  invariant 
with  respect  to  spatial  translation, 

(e(ri)e(r2))  =  C{ti  -  T2).  (3.17) 

One  may  argue  that  stationarity  can  never  be  achieved  in  a  real  environment.  Nev¬ 
ertheless,  from  a  practical  point  of  view,  it  is  a  reasonable  assumption  if  we  conflne 
ourselves  to  a  small  enough  region  so  that  the  change  across  the  whole  region  is  little. 
Alternatively,  a  random  process  can  be  characterized  in  the  wavenumber  domain  by 
its  power  spectrum 


/oo  roo  /“OO 

/  /  C'(r)  exp(— ik  •  r)dr,  (3.18) 

-OO  J—oo  J—00 

where  k  =  {kx,  ky,  kz)  is  the  wavenumber  vector. 

For  a  horizontal  correlation  function  to  be  Gaussian  or  exponential,  we  assume 
that  the  correlation  function  of  the  quantity  e(r)  is  horizontally  and  vertically  decou¬ 
pled  [7,  28,  35]: 


(^(ri  -  rz)  =  (e(ri)e(r2))  =  a^N{l)M{\zi  -  zz]),  (3.19) 

where  the  horizontal  vector  1  =  [(xi  —  X2),  {yi  —  2/2)],  and  a  is  the  standard  deviation 
of  €.  Geoacoustic  evidence  [29]  showed  that  the  vertical  correlation  length  Iz  is  much 
smaller  than  the  horizontal  correlation  length  lx  and  ly.  For  the  sake  of  convenience, 
the  vertical  correlation  function  M{z)  is  chosen  to  be  exponential,  i.e., 

M{z)  =  exp(-^).  (3.20) 
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Its  Fourier  transform  is 


The  horizontal  correlation  function  N{1)  can  be  both  Gaussian  and  exponential.  For 
a  Gaussian  correlation  function,  we  have 

N{x,y)  =  exp(-i[(^)2  +  (^)2]).  (3  22) 

Z  lx  ty 

(Notice  a  factor  of  1/2  difference  in  the  exponential  compared  to  that  in  Chernov’s 
[33]  and  Tang’s  [35]  work.)  The  corresponding  Fourier  transform  is 

p/2  I  p/2 

Wik,,  ky)  =  27rU,exp(-^-^  ^  ^).  (3.23) 

As  for  the  exponential  correlation  function,  there  are  two  choices;  one  is  for  an 
isotropic  random  field. 


Nil)  =  exp(-^), 
‘0 


(3.24) 


with  Zq  the  correlation  length  and  the  Fourier  transform 


W(k,)  = 


the  other  is  for  an  anisotropic  random  field 


N{x,y)  =  exp(-[J^  +  M]) 


Zx  Zy 


with  a  corresponding  Fourier  transform 


Wik,,ky)  = 


4ZxZy 


(3.25) 


(3.26) 


(3.27) 


ii  +  kiiDii  +  k^il)- 

The  second  function  will  not  degenerate  to  the  first  one  even  if  =  I  because 
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of  a  discontinuity  at  zero  lag.  However,  as  discussed  by  Ivakin  and  Lysanov  [28], 
this  can  be  corrected  by  introducing  a  second  length  scale  to  represent  the  minimum 
inhomogeneity  scale  in  Eq.3.24,  thereby  obtaining 

exp(— 1^)  —  eexp(— 1^) 

iV(l)  =  ^  (3-28) 

where  e  =  h/lo  and  li  Iq.  Under  the  condition  that  e  — >  0,  it  is  equivalent  to 
Eq.3.24. 

The  Von  Karman  function,  which  has  a  power-law  type  of  spectrum,  is  used  as 
a  correlation  function  to  describe  turbulence  (Tatarski  [65]),  to  characterize  hetero¬ 
geneity  in  the  earth’s  crust  (Frankel  and  Clayton  [63]),  and  to  model  the  seafloor 
morphology  (Goff  and  Jordan  [66]).  In  this  work,  we  will  follow  Goff  and  Jordan’s 
approach.  First,  a  set  of  functions  are  defined 

GM^))  =  b''ir)K^{bir))  0  <  6  <  oo,  (3.29) 

where  K,,  is  the  modified  Bessel  function  of  order  u  and 


5(r)  =  [v^QvY^'^.  (3.30) 

Here  (T)  stands  for  the  transpose  operation,  and  Q  is  a  positive-definite,  symmetric 
matrix  whose  Cartesian  elements  qij  have  dimensions  of  {length)~^.  In  the  three- 
dimensional  case,  it  can  be  expressed  in  terms  of  its  normalized  eigenvectors  On,  Cg, 
and  et  and  its  ordered  eigenvalues  k'^>k1>  kf, 

Q  =  kle^el  +  k'^CseJ  -k  k^ete^.  (3.31) 

The  Q  matrix  carries  information  about  the  anisotropy  of  the  random  field.  While 
kn,  kg  and  kt  correspond  to  the  inverses  of  the  correlation  lengths  in  three  dimensions 
respectively.  For  a  two-dimensional  case,  simply  drop  the  kt  term  in  the  expression. 
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The  correlation  function  can  then  be  set  to  be 


C{t) 


Mbjr)) 
G.(0)  • 


The  power  spectrum  for  a  two-dimensional  random  field  is 


(3.32) 


Wik:,,ky) 


A'Kva'^lxly 


(1  -b  kill  -b  kiiiY+^ 
A'KoHjy 

[  {I  +  kill  +  kill) 


(3.33) 


Because  we  are  interested  in  monostatic  backscattering  in  this  work,  azimuthally- 
summed  random  fields  will  be  generated  as  shown  in  Chapter  4.  The  x  and  y  axis  are 
then  chosen  to  be  overlaid  with  On  and  Cs,  i.e.,  no  skew.  Figures  3-1  and  3-2  show 
the  correlation  functions  and  power  spectra  for  several  value  of  i>  (same  as  Fig.  2  in 
Ref.  [66]).  An  asymptotic  roll-off  rate  -2{iy  +  1)  of  the  spectrum  is  related  to  the 
fractal  (Hausdorff)  dimension  D fractal  of  such  a  stochastic  process  as; 


D  fractal  —  3  U. 


(3.34) 


Therefore,  decreasing  the  parameter  u  would  increase  the  roughness,  which  is  the 
case  for  the  three-dimensional  case  as  well.  In  the  special  case  ol  u  =  1/2,  the 
power  spectrum  reduces  to  the  isotropic  exponential  power  spectrum  when  l^.  =  ly. 
Of  course,  we  can  construct  a  vertically  uncorrelated  random  field  model  with  the 
randomness  in  each  layer  specified  by  the  above  two-dimensional  power-law  type  of 
power  spectrum.  For  a  correlated  three-dimensional  case,  the  power  spectrum  is 

Wik  k  k  )  =  +  3/2) 

1  X,  y,  z)  ^2/2  +  f,2l2  +  jt2/2)^+3/2  >  (3-35) 

where  F  represents  the  Gamma  function.  The  roll-off  rate  now  is  -2(i^  -b  3/2)  and 
the  related  fractal  (Hausdorff)  dimension  is 

D  fractal  =  4  —  1^.  (3.36) 


86 


Figure  3-1:  Functional  form  of  the  correlation  function  G^ib)  plotted  for  values  of 
=  0, 1/2,  and  1. 


Figure  3-2:  Normalized  power  spectra  plotted  in  log-log  space. 
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Figure  3-3:  Normalized  isotropic  power  spectra  for  different  kinds  of  correlation  func¬ 
tions  (lx  —  ly  =  Im,  1/  =  1  in  power  law  distribution). 

As  can  be  seen  in  Fig.  3-3,  the  isotropic  power  spectra  (lx  =  ly)  for  Gaussian, 
exponential  and  power  law  distributions  have  different  roll-off  rates  versus  wavenum¬ 
ber  even  with  the  same  correlation  length.  For  a  Gaussian  distribution,  the  spectral 
level  at  high  wavenumber  is  low,  which  translates  to  a  lack  of  small-scale  roughness. 
On  the  other  hand,  there  are  significant  amounts  of  energy  in  high  wavenumbers  for 
exponential  and  power  law  distributions.  This  is  consistent  with  the  fact  that  there 
is  abundant  small-scale  roughness  with  the  latter  two  distributions. 


3.3  Backscattering  from  Sound  Speed  Variations 
in  the  Free-Space  Scenario:  Analytical  Solu¬ 
tions 

Analytical  solutions  for  volume  backscattering  can  be  achieved  in  some  special  situ¬ 
ations  such  as  the  free-space  scenario.  By  free  space,  we  mean  that  the  background 
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Free-Space  Scenario 


Figure  3-4:  The  free-space  monostatic  backscattering  scenario. 


sound  speed  and  density  are  constants  throughout  the  media.  Because  monostatic 
backscattering  is  the  concentration  of  this  work,  solutions  will  be  given  only  for  this 
scenario.  Meanwhile,  only  sound  speed  variations  will  be  discussed  for  the  sake  of 
mathematical  simplicity,  i.e.,  P  =  0(the  effect  of  density  variations  will  be  studied 
later  in  Chapter  5).  Figure  3-4  shows  the  scenario  that  we  are  going  to  study  in 
this  section.  In  the  figure,  H  represents  source/receiver  height,  6  is  the  scattering 
angle,  D  is  the  thickness  of  the  random  scattering  layer,  and  R  is  the  distance  be¬ 
tween  source/receiver  and  scatterers.  We  choose  the  background  sound  speed  to  be 
1500m/s  and  the  density  to  be  1000/c^/m®  for  the  rest  of  the  chapter. 


In  free  space,  the  Green’s  function  for  a  point  source  is  well  known  as: 


G(Rr,r)  =  - 


exp(— zk  •  R) 


(3.37) 


where  R  =  Rr  —  r  and  R  =  |R|.  Therefore,  for  monostatic  backscattering  and  unit 
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source  strength,  Eq.3.15  becomes 


If  we  have  Ri  and  R2  much  larger  than  the  scattering  volume  dimension  L,  with 
L  much  larger  than  the  correlation  length  of  the  random  field,  Ri  and  R2  can  be 
approximated  by  the  distance  between  the  source/receiver  and  a  reference  point  in 
the  scattering  volume  R.  After  performing  the  coordinate  transformation  to  center 
coordinate  Tq  and  relative  coordinate  we  have 


"  (feSF  I L  *»■  (3.39) 

Bear  in  mind  that  the  correlation  function  depends  only  on  the  relative  coordinate. 
By  expanding  the  integration  limit  of  the  relative  coordinate  to  infinity,  we  have 
actually  neglected  the  infiuence  due  to  the  existence  of  physical  boundaries,  which 
will  result  in  some  differences  when  compared  with  solutions  including  the  boundary 
effect.  For  a  thin  layer  of  scatterers,  the  equivalent  surface  backscattering  coefficient 
E  can  be  obtained  by 


E  = 


IsR'^ 

linA^ 


(3.40) 


where  is  the  intensity  of  the  incident  source  wave  a  unit  distance  away  from  the 
scattering  region,  R  is  defined  above  as  the  distance  between  the  source/receiver  and 
a  reference  point  in  the  scattering  region,  and  A  is  the  insonified  area(here  actually 
the  top  of  the  scattering  volume)  [52].  The  backscattering  strength  is  simply 


=  lOlogipE. 

For  a  point  source  with  unit  source  strength, 

/.  = _ ^ 

(47ri?)2- 


(3.41) 


(3.42) 
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Therefore,  E  can  be  calculated  by 


S  =  ^  C{i)  exp(-2ik  •  dro-  (3.43) 

In  the  rest  of  the  section,  we  will  compute  E  for  different  correlation  functions  that 
have  been  introduced  in  the  previous  section.  Although  the  goal  of  this  work  is  to 
model  monostatic  backscattering  from  the  inhomogeneities  of  a  full  three-dimensional 
random  field,  we  opt  to  investigate  first  backscattering  from  two-dimensional  pseudo¬ 
surfaces  with  sound  speed  variations  (instead  of  height  perturbations  in  the  usual 
surface  scattering  problem)  in  hopes  of  obtaining  some  insights  before  complicating 
the  problem  further  with  another  dimension.  In  fact,  scattering  from  this  pseudo¬ 
surface  is  equivalent  to  that  from  a  vertically  uncorrelated  Im  layer  with  the  same 
horizontal  correlation  structure.  While  at  the  end  of  this  section,  the  effect  of  vertical 
correlation  will  be  evaluated. 


3.3.1  Backscattering  from  2-D  Isotropic  Random  Field 

By  isotropic  we  mean  that  the  correlation  lengths  for  the  random  field  in  x  and  y 
direction  and  ly  are  equal  to  Ir- 

For  a  Gaussian  correlation  function,  the  backscattering  coefficient  can  be  obtained 
as: 

E(d)  =  — - — -  exp{—2k'^l^  cos^  9),  (3-44) 

ZTT 

where  9  is  the  scattering  grazing  angle. 

Figure  3-5  shows  the  relation  between  the  backscattering  strength  and  the  product 
of  horizontal  wavenumber  K  and  a,  which  equals  to  \/2lj.  for  the  Gaussian  distribu¬ 
tion.  It  can  be  noticed  that  the  backscattering  strength  has  a  peak  around  Ka  =  1 
and  decays  rapidly  otherwise.  As  discussed  by  Chernov  [33],  A'a  «C  1  corresponds 
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Figure  3-5:  Backscattering  from  sound  speed  variations  with  a  2-D  isotropic  Gaussian 
correlation  function. 

to  backscattering  from  small-scale  inhomogeneities,  which  results  in  isotropic  scat¬ 
tering.  For  Ka  »  1,  the  scattered  power  will  be  concentrated  at  a  small  forward 
angle;  therefore  the  backscattering  is  weakened.  The  reason  that  the  backscattering 
has  only  one  “resonance”  for  the  Gaussian  distribution  has  a  lot  to  do  with  the  fact 
that  the  Gaussian  distribution  is  a  single-scale  model. 

As  for  an  isotropic  exponential  correlation  function,  the  backscattering  coefficient 
is 


27r[l-h(2A;/,cos0)2]3/2  (3-45) 

As  shown  in  Fig.  3-6,  the  backscattering  strength  decays  much  slower  versus  Ka 
than  with  a  Gaussian  distribution.  Here  K  again  is  the  horizontal  wavenumber  and 
a  Ij..  The  explanation  for  this  behavior  is  that  a  random  field  with  an  exponential 
correlation  function  has  its  energy  distributed  almost  equally  among  all  wavenumbers, 
i.e.,  it  is  a  multi-scale  model.  It  means  that  the  random  field  would  have  roughness 
at  all  the  scales,  which  makes  it  an  effective  scatterer  for  various  wavelengths  or 
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Figure  3-6:  Backscattering  from  sound  speed  variations  with  a  2-D  isotropic  expo¬ 
nential  correlation  function. 

scattering  angles.  We  would  expect  that  the  backscattering  strength  has  a  different 
decaying  rate  versus  scattering  grazing  angle  for  a  Gaussian  correlation  function  and 
an  exponential  correlation  function  as  well. 


For  a  correlation  function  with  a  general  power-law  type  of  power  spectrum,  we 
have  the  following  backscattering  coefficient: 


m  = 


< 


V 


Trfl  -V  (2klr  cos 
7r[l  -f-  {2klr  cos  Oy] 


(3.46) 


Figure  3-7  depicts  that  the  backscattering  strength  versus  Ka  ioi  v  =  1  with 
a  =  Ir-  Similar  to  the  random  field  with  an  exponential  correlation  function,  the 
backscattering  strength  decays  slowly  when  Ka  increases. 


If  we  plot  together  the  backscattering  strength  for  the  above  three  types  of  corre¬ 
lation  functions,  it  is  evident  that  the  decaying  rate  versus  Ka  is  essentially  related 
to  the  roll-off  rate  of  their  power  spectra.  It  is  the  roll-off  rates  that  determines  the 
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Figure  3-7:  Backscattering  from  sound  speed  variations  with  a  2-D  isotropic  Von 
Karman  function  as  correlation  function  {u  =  1). 


Figure  3-8:  Backscattering  from  sound  speed  variations  with  different  kinds  of  2-D 
isotropic  correlation  functions. 
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roughness  at  small  length  scales  therefore  the  backscattering  at  large  Ka. 


3.3.2  Backscattering  from  2-D  Anisotropic  Random  Field 

The  effect  of  anisotropy  on  monostatic  backscattering  will  be  studied  in  the  following. 
By  anisotropic  we  mean  that  the  correlation  length  lx  ^  ly  Since  the  scenario  with  a 
co-located  point  source  and  receiver  is  of  interest  in  this  work,  an  annulus  would  be 
insonified  for  each  scattering  angle. 

In  the  case  of  a  Gaussian  correlation  function,  the  backscattering  coefficient  can 
be  derived  as 

1  /-Stt 

=  — A  2  ^  I  exp(— cos^  9{ll  cos^  4>  +  ll  sin^  <i>))d(f> 

Tt/T  •/  0 

i.4^2i  I 

=  cos^  0{ll  +  ll))Io{k^  cos^  0{ll  -  ID),  (3.47) 

where  <f)  is  the  azimuthal  angle  and  Iq  represents  the  zeroth  order  modified  Bessel 
function  of  the  first  kind.  Figure  3-9  shows  the  backscattering  strength  versus  graz¬ 
ing  angle  for  different  combinations  of  horizontal  correlation  lengths  in  the  x  and 
y  directions.  It  appears  that  slopes  of  backscattering  strength  versus  grazing  angle 
are  different  for  different  combinations.  In  the  isotropic  case,  a  random  field  with 
Ir  =  0.6m  is  more  effective  in  backscattering  than  that  with  Ir  =  Im  and  further  with 
Ir  =  2m.  Therefore,  on  the  one  hand,  a  random  field  with  ly  =  0.6m  would  improve 
the  overall  backscattering  ability  compared  to  one  with  ly  =  Im.  On  the  other  hand, 
a  random  field  with  ly  =  2m  would  have  a  negative  impact  on  total  backscattering. 
The  azimuthal  dependence  cannot  be  discerned  because  of  the  omnidirectional  source 
and  receiver. 

For  an  anisotropic  exponential  correlation  function,  the  backscattering  coefficient 
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Figure  3-9:  Backscattering  from  sound  speed  variations  with  a  2-D  anisotropic  Gaus¬ 
sian  correlation  function. 


is 


m  = 


r27r _ ^ _ 

27r2  Jo  [1  -f  (2klx cos 6 cos  4>Y][l-\-  {2kly  cos 9  sin  <j)Y] ' 


(3.48) 


What  we  see  in  Fig.  3-10  is  similar  to  that  in  Fig.  3-9.  However,  the  differences 
among  three  combinations  are  not  as  big  as  that  for  a  Gaussian  distribution.  From 
Fig.  3-8,  one  can  notice  that  the  decay  rate  versus  Ka  is  much  smaller  for  an  ex¬ 
ponential  distribution  than  for  a  Gaussian  distribution.  As  a  result,  the  change  of 
correlation  length  in  one  direction  will  not  have  as  strong  an  influence. 


For  an  anisotropic  correlation  function  specified  by  Eq.3.32,  the  backscattering 
coefficient  becomes 


E(0)  = 


r  k'^ua'^ljy  r27r 

27r2  Jo  [1  +  {2klj:  cos  9  cos  0)2  -f  (2kly  cos  9  sin  0)2]'^+i  ^  ^  ^ 


k^oHjy  /•2’r 


/ 


d(l> 


(3.49) 


i/  =  0 


27r2  Jo  [1  -1-  (2/c/a;  cos  9  cos  0)2  +  [2kly  cos  9  sin  0)2] 

As  can  be  seen  in  Fig.  3-11,  the  differences  in  backscattering  strength  due  to  the 
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Figure  3-10:  Backscattering  from  sound  speed  variations  with  a  2-D  anisotropic  ex¬ 
ponential  correlation  function. 


Figure  3-11:  Backscattering  from  sound  speed  variations  with  a  2-D  anisotropic  Von 
Karman  function  as  the  correlation  function  {1/  =  1). 
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change  of  correlation  length  in  the  y  direction  are  larger  than  that  with  an  exponential 
distribution  but  smaller  than  that  with  a  Gaussian  distribution.  This  is  consistent 
with  the  fact  that  the  decay  rate  of  the  backscattering  strength  versus  Ka  in  this  case 
is  larger  than  that  of  an  exponential  distribution  but  smaller  than  that  of  a  Gaussian 
distribution. 


3.3.3  Effect  of  the  Vertical  Correlation 

For  the  sake  of  simplicity,  the  random  sound  speed  and  density  variations  are  some¬ 
times  modeled  as  vertically  uncorrelated  [7].  What  we  will  show  next  is  the  difference 
in  backscattering  strength  versus  grazing  angle  for  vertically  correlated  and  uncorre¬ 
lated  random  fields.  As  we  discussed  in  the  previous  chapter,  the  scheme  to  character¬ 
ize  volume  backscattering  using  equivalent  surface  scattering  strength  versus  grazing 
angle  is  only  valid  when  the  scattering  layer  is  thin  with  respect  to  the  wavelength 
and  the  source/receiver  are  far  away.  The  environmental  parameters  will  be  chosen 
to  fulfill  this  requirement.  In  the  following,  we  will  present  only  analytic  solutions  for 
backscattering  from  random  fields  with  an  isotropic  Gaussian  or  exponential  distri¬ 
bution  horizontally  and  exponential  distribution  vertically. 


With  an  exponential  correlation  function  in  the  vertical  direction,  the  backscat¬ 
tering  coefficient  for  a  horizontally  isotropic  Gaussian  distribution  is 


^4^2/2/  rt 


(3.50) 


where  4  is  the  vertical  correlation  length  and  D  is  the  layer  thickness.  For  a  horizon¬ 
tally  isotropic  exponential  distribution,  the  backscattering  coefficient  becomes 


m  = 


k^aH%D 


7r[l  +  {2kl,  sin  0)2] [1  +  {2klr  cos  0)2]3/2  • 


(3.51) 


It  is  evident  from  the  above  two  equations  that  vertical  correlation  does  affect  the 
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Figure  3-12:  Backscattering  from  sound  speed  variations  with  an  isotropic  Gaussian 
distribution  horizontally  and  an  exponential  correlation  vertically  (cr  =  0.05,  Ix  —  ly  — 
2m,  /  =  400/^2,  H  —  400m,  D  =  7.5m). 

angular  dependence  of  the  backscattering  strength.  The  effect  is  insignificant  when 
{2klzmv9Y  C  1,  which  can  be  approximated  by  a  vertically  uncorrelated  model. 
However,  this  condition  can  be  satisfied  only  at  very  low  frequency  and  small  ver¬ 
tical  correlation  length.  Otherwise,  such  an  effect  should  not  be  ignored.  When 
{2klzSm6Y  1,  small  changes  in  Iz  will  not  have  strong  influences  on  the  angu¬ 
lar  dependence  but  only  the  overall  levels  of  backscattering  strengths.  While  with 
{2klzS\ndY  «  1,  we  would  expect  that  the  angular  dependence  of  the  backscattering 
strength  to  be  very  sensitive  to  the  change  of  klz- 

Figure  3-12  and  Fig.  3-13  are  two  examples  to  show  the  effect  of  vertical  correla¬ 
tion.  Figure  3-12  depicts  the  backscattering  strength  for  a  7.5m  thick  random  layer 
with  a  Gaussian  correlation  horizontally.  Figure  3-13  is  for  an  exponential  distribu¬ 
tion  horizontally.  The  correlation  lengths  4  and  ly  are  chosen  to  be  both  2m.  The 
standard  deviation  of  sound  speed  variations  is  5%  and  the  frequency  is  400Hz.  We 
present  the  backscattering  strength  curves  for  different  vertical  correlation  lengths 
Iz-  In  both  plots,  there  is  a  thin  solid  line  representing  the  uncorrelated  case,  for 
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Figure  3-13:  Backscattering  from  sound  speed  variations  with  an  isotropic  exponential 
distribution  horizontally  and  an  exponential  correlation  vertically  (cr  =  0.05,  l^  =  ly  = 
2m,  f  —  400Hz,  H  =  400m,  D  =  7.5m). 

which  we  use  the  solutions  of  a  2-D  isotropic  distribution.  The  justification  for  this 
comes  from  the  fact  that  the  3-D  solution  will  not  converge  to  the  uncorrelated  case 
with  1^  =  0  due  to  the  discontinuity  of  the  exponential  correlation  function  discussed 
earlier.  Fortunately,  what  interests  us  most  here  is  the  angular  dependence  of  the 
backscattering  strength,  which  is  well  embedded  in  the  2-D  solution.  As  can  be 
seen  from  both  figures,  the  angular  dependence  of  the  backscattering  strength  for 
Iz  =  0.07m  is  very  close  to  that  of  an  uncorrelated  model.  However,  as  increases, 
we  can  see  the  backscattering  strength  level  increases  for  smaller  grazing  angles,  which 
cannot  be  predicted  by  a  vertically  uncorrelated  model.  In  Fig.  3-8,  the  backscat¬ 
tering  strengths  possess  peak  values  for  certain  Ka.  Similar  behavior  is  expected  for 
the  backscattering  strength  versus  kzh,  where  kz  is  the  vertical  wavenumber.  The 
backscattering  strength  includes  contributions  from  both  the  horizontal  correlation 
and  vertical  correlation.  As  the  grazing  angle  decreases,  the  horizontal  wavenumber 
K  increases  and  the  vertical  wavenumber  kz  decreases.  The  contribution  due  to  the 
horizontal  correlation  will  decrease  since  Ka  increases,  i.e.,  moving  away  from  the 
peak,  while  the  contribution  from  the  vertical  correlation  might  increase  because  of  a 
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smaller  kj^.  As  a  result,  there  is  a  possibility  for  the  backscattering  strength  curve  to 
ascend  at  small  grazing  angles,  which  would  not  happen  for  a  vertically  uncorrelated 
random  medium. 

In  conclusion,  modeling  the  vertical  correlation  of  the  random  sound  speed  and 
density  variations  is  important  in  order  to  capture  the  full  scope  of  the  scattered  field, 
which  has  justified  our  efforts  to  model  backscattering  from  3-D  random  volume  in¬ 
homogeneities  in  the  sediment. 


101 


Chapter  4 


Generation  of 

Azimuthally-Summed  3-D  Random 
Sound  Speed  Variations 

4.1  Statement  of  the  Problem 

The  objective  of  this  research  is  to  develop  a  model  that  can  account  for  monostatic 
backscattering  from  3-D  volume  inhomogeneities. 

The  scenario  of  a  water  layer  on  top  of  the  sediment  layers  is  considered  here  with 
source/receiver  co-located  in  the  water  column.  By  applying  the  method  of  small 
perturbation  and  the  Born  approximation,  the  scattered  field  P®  due  to  sound  speed 
perturbation  =  0  here  for  simplicity)  in  the  sediment  volume  can  be  expressed  as 
(cf.  Eq.3.14): 


P\u,  R)  =  Am  [  e(Ri)(2  klG^R,  Rj))  dVi. 

where  A  is  the  sound  strength  and  m  is  the  density  ratio  between  the  top  of  the  sedi¬ 
ment  layer  and  the  water,  lo  is  the  radial  frequency,  R  and  Ri  are  the  source/receiver 
position  and  scatterer  position,  respectively.  The  relative  sound  speed  perturbation 
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is  represented  by  e  =  (5c/co,  ki  is  the  background  wavenumber  in  the  sediment,  Vi 
is  the  scattering  volume  and  G(R,Ri)  stands  for  the  Green’s  function  between  the 
scatterer  and  the  receiver.  If  the  above  integral  is  carried  out  in  cylindrical  coordi¬ 
nates,  dVi  can  be  expanded  as  r  dr  dcj)  dz.  By  observing  the  above  integral,  it  can  be 
seen  that  only  e  is  dependent  on  (j).  So  instead  of  generating  e  in  3-D,  if 

r2-jr 

r}{r,z)=  €{r,<j),z)d(p,  (4.1) 

J  0 

can  be  generated,  the  3-D  problem  is  reduced  to  a  2-D  problem.  Since  the  scattering 
intensity  is  of  most  interest  to  us,  the  generated  sound  speed  perturbation  should 
possess  the  desired  second-order  statistics. 

As  for  the  random  field  generation,  the  algorithms  can  be  catergorized  into  two 
classes:  one  that  is  implemented  in  the  spatial  domain  and  another  that  is  imple¬ 
mented  in  the  spectral  domain.  The  spectral  method  [19]  is  more  popular  than 
the  spatial  method  in  the  underwater  acoustics  community  because  of  its  efficiency 
in  generating  large,  multi-dimensional  fields.  It  was  applied  by  Thorsos  [67]  in  his 
study  of  rough  surface  scattering  from  a  1-D  Gaussian  random  surface.  The  2-D 
Goff-Jordan  power  law  spectrum  [66]  is  also  extensively  used  to  simulate  the  seafloor 
[7,  68].  Nevertheless,  the  literature  on  3-D  random  field  generation  can  hardly  be 
found  in  underwater  acoustics,  in  contrast  to  the  large  number  of  papers  in  water 
resource  research  literature  [69,  70,  71].  The  reason  may  rest  on  the  fact  that  the 
computational  requirement  is  too  big  for  the  scale  of  the  random  field  that  we  are 
dealing  with  in  underwater  acoustics.  To  serve  the  purpose  of  this  research,  the  spec¬ 
tral  method  is  used  because  of  its  simplicity  and  the  availability  of  analytical  forms  of 
different  power  spectral  density  functions  to  describe  the  sound  speed  perturbation. 
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4.2  Formulas 


The  spectral  method  for  generating  an  azimuthally-summed  3-D  random  field  r]{r^  z) 
will  be  presented  in  this  section.  For  the  sake  of  simplicity,  generation  of  azimuthally- 
summed  2-D  random  fields  will  be  studied  first.  In  other  words,  the  dimension  z 
(corresponding  to  the  vertical  direction)  is  dropped  for  the  time  being.  However,  it 
can  be  easily  added  to  the  formulas  later. 


In  Cartesian  coordinates,  the  power  spectral  density  W  can  be  obtained  by 


W{k^,ky)  =  lim 

L-^oo  ij^ 


U 


Ll‘1  /•L/2 


L/2  J-L/2 


2\ 

^‘^^dxdy  V 


(4.2) 


where  k^,  ky  are  the  wavenumbers  in  the  re,  ^  directions,  is  the  patch  size  assuming 
square  area,  ( • )  stands  for  the  ensemble  average  operation  and  |  •  |  indicates  the 
modulus  operation.  The  power  spectral  density  should  be  normalized  such  that 


roo  foo 

(2^)2  J  J  ^^{kxt  ky)  dkx  dky  =  g  ,  (4-3) 

where  o  is  the  standard  deviation  of  the  sound  speed  perturbation.  Equation  4.2  can 
be  proved  to  satisfy  the  relation  [72] 


/oo  roo 

/  (4.4) 

-OO  ^—OO 

where  C{ix,  iy)  is  the  correlation  function  and  ^y  are  the  lag  distances  between 
two  points  in  the  x  and  y  directions.  The  consistency  between  Eq.4.2  and  Eq.4.4 
has  an  important  implication.  It  guarantees  that  if  the  power  spectral  density  of 
the  generated  random  field,  estimated  through  Eq.4.2,  converges  to  the  desired  one, 
the  generated  random  field  would  have  the  desired  correlation  function  according  to 
Eq.4.4,  and  therefore  the  desired  second-order  statistics. 
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According  to  Eq.4.2,  a  random  field  can  be  generated  by 


<^.y)  =  r  r  x/WK(A:x,  ky)L^Rn{k.,  dky.  (4.5) 

J-oo  J-OO  '' 

Here  the  random  number  Rn  fits  the  distribution  Rn{kx,  ky)  ~  [A^(0,  l)+i  A^(0,  l)]/v^, 
where  A^(0, 1)  stands  for  a  Normal  distribution  with  zero  mean  and  standard  deviation 
1.  In  order  for  t{x,y)  to  be  real,  Rn{kx,ky)  has  to  be  conjugate  symmetric  between 
the  first  and  the  third  quadrants  as  well  as  the  second  and  the  fourth  quadrants.  If 
the  following  coordinate  transformation 

{kx  =  kr  cos  6  \  X  =  r  cos  6 

and  I  ^  ,  (4.6) 

ky  =  krSm.6  I  y  =  r sin (j) 

is  performed,  Eq.4.5  will  become 

f(r,  4')  =  ^  j(°°  1^’  <JW{,k„0)V‘K{K,  dK  de.  (4.7) 

The  equation  for  r}{r)  will  be 

/•2jr  1  TOO  rltt  ! - 

yW{kr,6)L^Rn{kr,6)krJo{krr)dkrd9.  (4.8) 

We  point  out  that  rj{r)  is  a  nonstationary  process.  Here  A0  =  27r/N{r)  and  N(r)  is 
the  number  of  scatterers  around  the  circle  with  radius  r.  Assuming  constant  density 
for  the  scatterers,  N{r)  will  increase  proportionally  to  r.  Therefore,  the  variance  of 
r]{r)  will  decrease  accordingly,  which  makes  T]{r)  a  nonstationary  process.  Collecting 
the  items  on  the  right-hand  side  which  depend  on  0,  it  gives 

F{kr)  =  A  yJm^^^Rn{kr,  9)  d9 .  (4.9) 

Changing  the  continuous  integral  to  a  discrete  sum  and  using  the  relation 

=  Akx  Aky  =  kr  Akr  A9  {kr  +  0), 
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where  A(-)  represents  the  sampling  interval,  Eq.4.9  becomes 


M{kr) 

m)  =  E  ■ 

i=l 


Iwjkr,  9i)A9 
kr  Akj. 


Rn{kr,0i) 


(4.10) 


where  6i  =  iA6  and  M{kr)  is  the  number  of  samples  around  a  circle  of  radius  kr. 
The  value  of  F  at  kr  =  0  can  be  arbitrary  as  long  as  it  is  finite  because  it  will  be 
multiplied  by  kr,  which  is  0  here,  in  Eq.4.8.  F{kr)  can  be  treated  as  a  new  random 
variable  with  mean  0  and  variance 


M(kr) 

E 


i=l 


Wikr,0i)A0 

kr  Akr 


(kr  /  0), 


since  Rn{kr,  0i)  are  independent  of  each  other  and  each  has  the  variance  of  1.  There¬ 
fore,  F{kr)  can  be  replaced  by 


S’  w{k„e)de 

kr  dkr 


Rn{kr).  '{kr^Qi), 


where  Rn{kr)  ~  N{Q,  1).  Notice  that  here  Rn{kr)  is  a  real  number  due  to  the  fact  that 
the  imaginary  parts  of  Rn{kr,  0)  have  been  cancelled  because  of  the  earlier  mentioned 
property  of  conjugate  symmetry.  Eq.4.8  now  becomes 


.M  =  E^ 


W{kr„0)d0 
kr  Akr 


Rni,kri')kr^Jo{,krfT  'jAkr- 


(4.11) 


where  kn  =  i  Akr- 


As  mentioned  above,  the  right  second-order  statistics  of  generated  r]{r)  is  guar¬ 
anteed.  From  Eq.4.4,  the  correlation  function  of  e  in  Cartesian  coordinates  can  be 
written  as 

{^{xi,yi)e*{x2,y2))  = 

1  roo  roo 

y_oo  J-oc  dk,  dky.  (4.12) 
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where  {*)  stands  for  conjugate  operation.  If  the  coordinate  systerm  is  again  changed 
to  polar  coordinates,  it  will  result  in 


(e(n,0i)e*(r2,</>2))  = 

1  foo  /*27r 

-7^  /  (4.13) 

Jq  Jo 

Therefore,  the  correlation  function  of  r]{r)  is 

f2n  /•2n 

{v{n)v*{r2))  =  /  /  (e(ri,«^i)c*(r2,(^2))#i#2 

Jo  Jo 
foo  r  l■2■K 

=  U  W{kr,e)d6  krJo{krri)Jo{krr2)dkr.  (4.14) 
It  can  be  seen  that  r]{r)  generated  by  Eq.4.11  satisfies  the  above  equation.  Let  us 


denote 


f2Ti 

Ws=  Wikr,e)d9, 
Jo 


which  can  be  considered  the  radial  power  spectral  density.  It  is  instructive  to  compute 
the  spectral  estimates  Ws  from  the  random  realizations  of  sound  speed  perturbations 
for  comparison  with  Ws-  In  continuum  notation,  we  calculate 


/  fL/2  2, 

(  w{r)T]{r)rJo{krr)dr  j, 


(4.15) 


where  w{r)  is  a  real  non-negative  weighting  or  window  function  normalized  such  that 


O  rL/2 

—  /  w^{r)  dr  =  1. 

Lj  Jo 


(4.16) 


The  weighting  function  is  chosen  to  improve  the  spectral  estimation  by  reducing  an 
effect  known  as  spectral  leakage.  So  if  the  estimated  W  agrees  with  the  desired  Wg, 
the  generated  T]{r)  can  be  considered  satisfactory.  This  equation  will  be  used  later  in 
numerical  experiments  as  the  criterion  to  determine  the  acceptability  of  the  generated 
random  field. 
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Eq.4.11  can  be  easily  extended  to  the  generation  of  a  3-D  azimuthally-summed 
random  field  in  cylindrical  coordinates: 

1  oo  oo 
m=0n=-oo 

(4.17) 

where  =  mAkr,kz^  =  nAA:^  and  Rn{kr,k^)  has  to  be  conjugate  symmetric  in 
the  2:  direction  in  order  for  7y(r,  z)  to  be  real.  The  layer  thickness  is  denoted  by  D. 
This  process  can  be  viewed  as  passing  the  random  series  7;(r)  generated  according 
to  Eq.4.11  through  a  filter  in  the  z  direction.  As  for  the  power  spectral  density 
W{kr.,kz,0)  used  in  the  generation,  it  can  be  in  any  form,  including  the  result  of  a 
direct  measurement  in  an  experiment. 


If  the  area  within  a  few  wavelengths  distance  from  the  center  is  not  a  concern, 
Eq.4.17  can  be  evaluated  using  a  Fourier  transform  instead  of  a  Hankel  transform  in 
the  r  direction  and  a  Fourier  transform  in  z  direction,  which  is  essentially  the  same 
as  the  Fast-Field-Program(FFP)  approach  introduced  by  DiNapoli  and  Deavenport 
[73].  First  the  Bessel  function  Joikrv)  is  expressed  in  terms  of  Hankel  functions, 

Mkrr)  =  i  ,  (4.18) 

where  and  are  the  zeroth  order  Hankel  functions  of  the  first  and  second 
kind,  respectively.  Next,  HQ^\krr)  and  H^\krT)  are  replaced  by  their  asymptotic 
forms  [74] 


lim  H^\krr)  = 
lim  H^\krr)  = 

fcrr->oo  “  ^ 


(4.19) 
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to  arrive  at  the  following  expression, 


Ah  I  Ah  °°  °°  r  I  r27r 

.  (4.20) 


In  the  model/data  comparison  presented  in  Chapter  6,  our  interest  is  on  medium  to 
low  grazing  angle  backscattering.  The  generated  random  field  close  to  the  center  is 
going  to  be  discarded  anyway.  As  a  result,  Eq.4.20  can  be  exploited  to  speed  up  the 
computation. 


4.3  Numerical  Examples 

In  this  section,  azimuthally-summed  2-D  random  fields  are  going  to  be  generated. 
These  numerical  examples  will  help  verify  the  methodology  that  is  described  above. 
The  reasons  that  2-D  random  fields  instead  of  3-D  random  fields  are  chosen  are  due  to 
the  following  considerations:  first  is  simplicity,  since  it  is  easier  to  show  the  agreement 
of  power  spectral  densities  and  correlation  functions  in  2-D  cases;  second  is  the  fact 
that  the  azimuthal  summation  is  performed  only  in  a  horizontal  plane. 

The  situations  considered  will  include  isotropic  and  anisotropic  random  fields.  By 
isotropic,  it  means  that  the  correlation  lengths  in  both  the  x  and  y  directions  are  the 
same;  by  anisotropic,  it  means  that  the  correlation  lengths  are  different  in  the  x  and 
y  directions.  In  both  cases,  Gaussian  correlation  functions  will  be  used.  However,  it 
can  be  regarded  as  a  general  assessment  of  the  methodology. 

The  2-D  Gaussian  correlation  function  in  Cartesian  coordinates  is 

C&,«,)=<7^exp(-i(|  +  |]),  (4.21) 

X  y 

where  a  is  the  standard  deviation,  and  lx  and  ly  are  the  correlation  lengths  in  the  x 


109 


and  y  directions,  respectively.  The  corresponding  power  spectral  density  is 


W{K,ky) 


o  /* 


2'KaHxly  exp(- 


(4.22) 


4.3.1  Isotropic  Case 

In  the  isotropic  case,  lx  =  ly  =  I-  By  performing  the  coordinate  transformation  from 
Cartesian  coordinates  to  polar  coordinates  according  to  Eq.A.3,  Eq.A.2  becomes 

P/2 

W  {kr,  6)  =  27ro-^/^  exp( — |-).  (4.23) 

Since  W  is  independent  of  0,  Wa{kr)  is  simply  27rW{kr).  The  correlation  function  of 
T){r)  can  be  obtained  as  (see  Appendix  A) 

{vin)v*{r2))  =  47rV/o(^)  exp(— (4.24) 

where  Iq  stands  for  the  zeroth  order  modified  Bessel  function  of  the  first  kind.  The 
variance  at  each  range  r  is 

ivWir))  =  47rV/o(^)  exp(-^).  (4.25) 

It  is  evident  that  the  variance  of  t]  depends  on  r,  which  proves  that  r){r)  is  a  nonsta¬ 
tionary  process. 

Figure  4-1  shows  10  realizations  of  the  generated  random  sound  speed  perturba¬ 
tions  after  being  summed  up  azimuthally.  Here  the  horizontal  correlation  lengths  lx 
and  ly  are  both  5m.  The  standard  deviation  cr  is  chosen  to  be  1%.  There  is  a  decaying 
trend  of  the  magnitude  over  range,  which  is  the  result  of  decaying  variances. 
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Figure  4-1:  10  realizations  of  azimuthally-summed  2-D  random  fields  with  a  Gaussian 
correlation  function  and  a  —  1%,  lx  =  ly  5m 


Figure  4-2:  Comparison  of  radial  power  spectral  density  estimate  with  exact  radial 
power  spectral  density  for  an  ensemble  of  100  realizations  (isotropic  case) 
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Figure  4-3:  Comparison  of  variance  estimate  with  exact  variance  for  an  ensemble  of 
100  realizations(isotropic  case) 

A  comparison  of  radial  power  spectral  density  estimate  Ws  and  the  exact  Ws  is 
shown  in  Fig.  4-2.  Here,  Wg  is  obtained  with  100  realizations  using  the  same  param¬ 
eters  as  in  Fig.4-1.  A  Hanning  window  is  chosen  to  be  the  weighting  function  w{r). 
Reasonable  convergence  of  Ws  to  Wg  is  observed  except  at  low  spectral  levels.  This 
difference  occurs  because  w{r)  still  has  a  small  finite  value  at  the  endpoints,  which  is 
a  typical  spectral  leakage  effect. 

Fig.  4-3  shows  a  comparison  of  the  variance  estimate  of  77  (r)  and  the  exact  result 
obtained  through  Eq.A.6.  The  estimate  is  obtained  by 

■J  ^rlz 

Var[n{r)]  =  — -  _  £  (4.26) 

where  Nriz  is  the  number  of  realizations  used  in  the  estimation.  Again,  100  realization 
has  been  used.  The  agreement  between  the  estimate  and  the  exact  result  is  reason¬ 
ably  good,  which  shows  the  feasibility  of  this  method.  More  ensemble  averaging  will 
certainly  improve  the  match. 
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Results  Obtained  by  the  Far-Reld  Approximation 


Figure  4-4:  Comparison  of  radial  power  spectral  density  estimate  with  exact  power 
spectral  density  for  an  ensemble  of  100  realizations  using  the  FFP  approach 

Fig.  4-4  is  the  same  as  Fig.  4-2  except  the  random  fields  are  generated  by  the 
FFP  approach  described  in  Eq.4.20.  The  good  fit  between  the  estimated  and  exact 
VFs  indicates  the  acceptability  of  this  approach. 


4.3.2  Anisotropic  Case 

In  an  anisotropic  situation,  the  coordinate  transformation  of  Eq.A.2  through  Eq.A.3 
will  give 


W{kr,6)  =  2na^lxly  exp{- 


cos^  ^  sin^  6) 


)■ 


(4.27) 


As  a  result, 


Ws{kr)  -  47rV^/o( - - — —)e 


(4.28) 
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Figure  4-5:  10  realizations  of  azimuthally-summed  2-D  random  fields  with  a  Gaussian 
correlation  function  and  a  =  1%,  =  2m,  ly  =  5m 


(See  Appendix  A  for  a  detailed  derivation).  The  variance  of  'q{r)  is 


Var{r]{r)'\  = 

r2w  r2‘K  ^2 

^  Jo  Jo  (4.29) 

X  y 

This  equation  can  be  evaluated  numerically. 


Figure  4-5  shows  10  realizations  of  the  generated  random  sound  speed  perturba¬ 
tions  after  being  summed  up  azimuthally.  Here  the  horizontal  correlation  lengths 
and  ly  are  2m  and  5m,  respectively.  The  standard  deviation  a  is  again  1%.  By  com¬ 
paring  Fig.4-5  with  Fig.  4-1,  one  can  notice  that  there  are  more  oscillations  over  the 
range  in  the  previous  one.  This  is  because  the  correlation  length  in  this  anisotropic 
case  is  2m  instead  of  the  5m  used  in  the  isotropic  case,  which  means  that  the  corre¬ 
lation  of  the  original  2-D  field  before  azimuthal  summation  is  not  as  good  as  that  in 
the  isotropic  case.  This  is  reflected  in  the  azimuthally-summed  field.  It  will  be  more 
evident  when  the  radial  power  spectral  densities  are  compared  later. 
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Figure  4-6:  Comparison  of  radial  power  spectral  density  estimate  with  exact  radial 
power  spectral  density  for  ensemble  of  100  realizations  (anisotropic  case) 

Again,  good  agreement  can  be  seen  between  the  estimated  power  spectral  density 
Wg  and  the  exact  Wg.  Wg  is  obtained  with  100  realizations  using  the  same  parameters 
as  in  Fig.4-5.  A  Hanning  window  is  chosen  to  be  the  weighting  function  w{r).  In  this 
particular  case,  the  difference  at  low  spectral  level  due  to  spectral  leakage  cannot  be 
observed  since  it  is  out  of  the  plotted  wavenumber  region.  The  match  between  the 
estimated  variance  of  T]{r)  and  the  exact  one  obtained  through  numerical  evaluation 
of  Eq.  A. 8  is  good  as  well. 

Next  we  compare  the  radial  power  spectral  density  Wg  for  random  fields  with 
different  correlation  lengths  lx-  It  can  be  seen  in  Fig.  4-8  that  the  decay  rate  of  power 
spectral  density  versus  wavenumber  is  inversely  proportional  to  the  correlation  length 
lx  when  ly  is  equal.  The  larger  the  lx,  the  better  the  correlation  of  the  original  field,  the 
larger  the  decay  rate.  From  the  perspective  of  spatial-wavenumber  relationship,  there 
will  be  more  energy  at  high  wavenumber  for  poorer  correlation,  which  corresponds  to 
higher  oscillation  in  spatial  domain  as  can  be  seen  in  Fig.  4-5  compared  to  that  in 
Fig.  4-1. 
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Figure  4-7:  Comparison  of  variance  estimate  with  exact  variance  for  an  ensemble  of 
100  realizations  (anisotropic  case) 


Figure  4-8:  Comparison  of  radial  power  spectral  densities  for  random  fields  with  the 
same  correlation  length  ly  but  different  correlation  length  The  estimated  radial 
spectral  densities  are  obtained  by  an  ensemble  average  of  100  realizations 
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4.4  Summary 


The  purpose  of  generating  azimuthally-summed  random  fields  is  to  reduce  the  compu¬ 
tational  load  of  studying  monostatic  backscattering  from  3-D  random  inhomogeneities 
by  transforming  a  3-D  problem  to  a  2-D  problem  if  the  source  and  receiver  are  co¬ 
located.  The  numerical  examples  given  in  the  previous  section  prove  that  the  proposed 
approach  is  successful,  which  will  facilitate  the  synthesis  of  backscattering  time  series 
in  the  next  chapter. 
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Chapter  5 


Numerical  Simulations  of 
Backscattered  Time  Series 


One  of  the  objectives  of  this  thesis  work  is  to  enable  the  simulation  of  the  backscat¬ 
tered  field  and  backscattering  time  series  due  to  three-dimensional  volume  inhomo¬ 
geneities  in  the  sediment.  In  Chapter  3,  a  small  perturbation  approach  was  employed 
to  derive  the  formulations  for  bottom  scattering  from  sound  speed  and  density  varia¬ 
tions.  In  Chapter  4,  a  scheme  was  proposed  to  generate  azimuthally-summed  three- 
dimensional  random  fields  using  the  spectral  method.  With  the  calculation  of  the 
Green’s  function  discussed  in  Chapter  2  handled  by  OASES  and  sound  speed  and 
density  variations  generated  by  the  above-mentioned  spectral  method,  we  are  pre¬ 
pared  to  simulate  monostatic  backscattered  fields  and  backscattering  time  series  at 
the  receiver. 

In  this  chapter,  we  will  describe  procedures  for  generating  backscattering  time 
series  first.  Then  the  emphasis  will  be  put  on  comparisons  between  numerical  simu¬ 
lation  and  analytic  solutions.  The  effect  of  density  variations  on  backscattering  will 
be  studied  at  the  end. 
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5.1  Numerical  Simulations 


The  keys  to  the  numerical  simulation  of  scattered  fields  and  the  scattering  time  series 
rest  on  the  following  equations  (cf.  Eq.3.14) 

P.(Rr)  =  I  e(r')[(2  +  /?)fcJ(/)Po(r')G(R„r')  -  /?(yPo(r')  •  VG(R„  r'))]*', 

where  Rr  stands  for  the  receiver  position  and  r'  stands  for  the  scatterer  position.  The 
scattering  time  series  y{t)  can  be  obtained  by  (cf.  Eq.3.16) 

1 

y{t)  =  —  Ps{Iir,u)S{u})exp{iujt)dh}, 

ZTT  J^oo 

where  u  is  the  radial  frequency  and  S(cv)  is  the  spectrum  of  the  source  signal.  The 
numerical  procedures  for  these  are  self-evident,  i.e., 

•  Calculate  the  Green’s  function  G  and  the  mean  field  Pq.  For  simple  cases 
such  as  a  free-space  scenario,  analytic  solutions  are  available,  while  for  more 
complicated  environments,  some  numerical  toolboxes  {e.g.,  OASES  and  Kraken 
[75])  can  be  exploited.  OASES  is  our  choice  in  this  work. 

•  Generate  azimuthally-summed  3-D  random  sound  speed  variations  using  the 
spectral  method  described  in  Chapter  4.  For  monostatic  backscattering,  ex¬ 
cept  for  e,  all  the  terms  on  the  RHS  of  Eq.3.14  are  independent  of  azimuth 
in  a  cylindrical  coordinate.  As  a  result,  we  can  take  advantage  of  having  an 
azimuthally-summed  random  field  instead  of  a  full  3-D  random  field  because  of 
the  difficulty  of  meeting  the  computational  requirement  for  the  generation  of 
the  latter. 

•  Obtain  the  monostatic  backscattered  field  using  Eq.3.14.  If  the  ultimate  goal 
is  to  simulate  a  backscattering  time  series,  scattered  fields  for  many  frequencies 
have  to  be  calculated.  Subsequently,  G  and  Pq  have  to  be  computed  for  mul¬ 
tiple  frequencies.  Meanwhile,  random  sound  speed  variations  need  only  to  be 
generated  once  because  they  are  frequency-independent. 
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Figure  5-1:  The  source  signal  and  its  spectrum. 

•  Produce  backscattering  time  series  by  Fourier  synthesis  of  the  scattered  field  in 
the  full  frequency  spectrum. 

Next  we  present  an  example  of  a  simulated  backscattering  time  series  at  one 
receiver.  We  will  utilize  the  same  free-space  scenario  as  shown  in  Fig.  3-4.  The  back¬ 
ground  sound  speed  is  chosen  to  be  1500m/s,  the  density  the  source/receiver 

height  400m  and  there  is  no  attenuation.  Figure  5-1  depicts  the  source  signal  that 
is  going  to  be  used  in  the  simulation  throughout  this  chapter  and  its  spectrum.  It  is 
a  Hanning-windowed  CW  wave  with  a  pulse  duration  of  0.1s.  The  sampling 

frequency  is  20{)QHz. 

Assuming  exponential  correlation  functions  horizontally  and  vertically  for  a  7.5m 
(2  wavelength)  thick  random  layer  with  the  horizontal  correlation  length  =  2m 
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Figure  5-2:  Simulated  backscattered  returns  from  an  inhomogeneous  layer.  The  cor¬ 
relation  function  of  random  sound  speed  variations  within  the  layer  are  exponential 
both  horizontally  and  vertically. 

and  vertical  correlation  length  1 2  =  0.5m,  we  can  simulate  backscattered  returns 
recorded  at  the  receiver  as  shown  in  Fig.  5-2.  The  single  source/single  receiver  config¬ 
uration  is  exploited  and  it  is  the  setup  for  the  rest  of  the  chapter.  The  source/receiver 
height  is  fixed  at  400m  from  now  on.  Only  sound  speed  variations  will  be  modeled 
for  the  time  being  and  the  effect  of  density  fluctuations  will  be  examined  in  the  last 
section  of  this  chapter.  Since  the  concentration  of  this  thesis  work  is  on  mid-  to 
low-  grazing  angle  backscattering,  we  have  tapered  the  random  field  close  to  normal 
incidence  with  half  of  a  Hanning  window.  This  causes  an  increasing  amplitude  at 
the  beginning  of  the  time  series.  The  amplitude  is  then  decaying  with  travel  time 
increasing  because  of  the  spreading  loss. 

Simply  by  observing  the  simulated  time  series,  it  is  almost  impossible  to  tell  if  the 
numerical  result  is  valid.  Therefore,  we  have  to  find  a  way  to  prove  that  the  generated 
backscattered  returns  make  sense,  which  will  be  the  task  of  the  next  section. 
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5.2  Comparisons  with  Analytic  Solutions  in  the 
Pree-Space  Scenario 

With  the  quantum  leap  in  the  computational  ability  in  recent  years,  numerical  for¬ 
ward  modeling  has  gained  more  and  more  popularity.  However,  there  is  always  a 
question  for  this  type  of  approach:  how  do  we  verify  the  numerical  results?  Usually 
there  are  two  ways  to  check:  the  first  and  ideal  scenario  is  to  make  comparisons  of 
numerical  results  with  appropriate  analytic  solutions,  pending  the  availability  of  such 
analytic  solutions;  the  second  is  to  check  if  the  model  under  development  agrees  with 
other  numerical  models  or  benchmark  solutions,  which  is  a  legitimate  verification  ex¬ 
cept  for  some  extreme  cases.  In  the  case  that  both  of  the  above  approaches  cannot 
be  carried  out,  the  last  but  less  convincing  approach  is  to  see  if  the  numerical  results 
make  sense  physically  in  solving  some  well-understood  problems.  Fortunately  for  this 
work,  we  have  analytic  solutions  for  monostatic  backscattering  from  volumetric  het¬ 
erogeneities  in  the  free-space  case,  which  makes  our  life  much  easier.  Formulations 
for  such  analytic  solutions  are  given  in  Chapter  3.  Subsequently,  analytic  results  will 
be  presented  in  this  chapter  without  being  reiterated  in  detail. 

As  we  know  in  numerical  modeling,  one  common  practice  is  to  discretize  the  con¬ 
tinuous  field,  i.e.,  sampling.  The  final  results  would  normally  vary  with  the  sampling 
rate  unless  adequate  sampling  of  the  field  has  been  achieved.  Therefore,  it  is  criti¬ 
cal  to  find  out  the  sufficient  sampling  rate,  thereby  reaching  a  unique  solution.  An 
effective  and  widely  used  procedure  is  the  convergence  test.  The  idea  is  to  increase 
the  sampling  rate  until  the  end  result  approaches  a  constant  value,*,  e.,  converges. 
The  lowest  sampling  rate  for  the  result  to  reach  the  constant  value  is  designated 
the  sufficient  sampling  rate.  Theoretically,  convergence  tests  need  to  be  performed 
for  every  problem.  Figure  5-3  shows  an  example  of  a  convergence  test.  The  ensem¬ 
ble  averages  of  backscattering  intensity  versus  time  are  plotted  for  50  realizations. 
Backscattering  intensity  here  is  defined  as  the  square  of  the  envelope  of  the  backscat- 
tered  returns.  Meanwhile,  from  a  statistics  point  of  view,  a  single  realization  is  not 
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Figure  5-3:  The  convergence  test  for  backscattering  from  2-D  sound  speed  variations 
with  a  power  law  distribution  (4  =  ly  =  Im). 

sufficient  to  determine  a  statistical  process.  Consequently,  ensemble  averaging  will 
instead  be  studied  for  the  rest  of  the  chapter.  The  result  shown  in  Fig.  5-3  is  for 
backscattering  from  2-D  sound  speed  variations  with  a  power  law  distribution.  The 
horizontal  correlation  lengths  lx  and  ly  are  Im,  respectively.  The  standard  deviation 
of  sound  speed  variations  a  is  5%.  From  the  figure,  it  can  be  seen  that  sampling 
interval  of  1/8  wavelength,  i.e.,  a  sampling  rate  of  8k  where  k  is  the  wavenumber, 
is  sufficient  for  the  backscattering  intensity  to  converge.  Of  course,  the  sufficient 
sampling  rate  will  largely  depend  on  what  type  of  correlation  function  we  have  for 
the  random  field,  or  more  specifically,  the  roll-off  rate  of  the  power  spectrum  ver¬ 
sus  wavenumber.  Convergence  tests  have  been  done  for  Gaussian  and  exponential 
distributions  as  well,  and  the  sampling  interval  of  1/8  wavelength  horizontally  has 
been  found  to  be  adequate,  although  results  are  not  shown  here  for  the  sake  of  brevity. 

As  mentioned  in  Chapters  2  and  3,  for  surface  scattering  or  volume  scattering  from 
a  thin  layer,  the  quantity  of  scattering  strength  versus  grazing  angle  can  be  used  to 
characterize  the  scattering  process.  There  is  no  angular  ambiguity  in  this  case  even 
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for  a  single  source/single  receiver  configuration  because  the  arrival  time  can  be  easily 
related  to  the  scattering  angle.  The  equation  for  the  calculation  of  the  backscattering 
strength  is  as  follows 


BSS  =101og,„^,  (5.1) 

where  /«  is  the  intensity  of  backscattered  signals  at  the  receiver,  is  the  intensity 
of  the  incident  sound  wave  a  unit  distance  away  from  the  scattering  area,  R  is  the 
distance  from  the  source  to  the  scatterer,  and  A  is  the  insonified  area  [52].  Following 
the  same  procedure  as  in  Ref.  [10],  we  have 

A  =  ncrR,  (5.2) 

where  c  is  the  sound  speed  and  r  is  the  pulse  duration.  At  the  same  time.  Is  and 
can  be  estimated  respectively  by 

4  =  ;  /  y\T)  dT,  (5.3) 

where  y{T)  is  the  scattered  return  at  the  receiver  and  t  is  the  round  trip  travel  time 
from  the  source/receiver  to  the  scatterers  and 

/in  =  -  /  s^{t)  dt,  (5.4) 

T  J  to 

where  S(t)  is  the  source  signal.  We  will  estimate  the  backscattering  strength  from 
the  simulated  backscattered  returns  and  compare  the  results  with  analytic  solutions 
obtained  in  Chapter  3. 

We  will  categorize  the  comparisons  into  the  2-D  isotropic  case,  the  2-D  anisotropic 
case  and  the  3-D  vertically  correlated  case. 
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Figure  5-4:  Backscattering  from  2-D  isotropic  sound  speed  variations  with  a  Gaussian 
correlation  function,  averaging  over  50  realizations  (/  =  400/^2;,  H  =  400m,  lx  =  ly  = 
Im,  a  =  0.05) 

5.2.1  2-D  Isotropic  Case 

First  we  will  present  an  example  of  backscattering  from  random  sound  speed  varia¬ 
tions  possessing  a  2-D  isotropic  Gaussian  correlation  function.  Figure  5-4  shows  the 
ensemble  average  over  50  realizations  of  simulated  backscattering  intensities  at  one 
receiver.  The  source/receiver  height  is  400m.  The  frequency  is  AOOHz.  The  horizon¬ 
tal  correlation  lengths  lx  and  ly  are  chosen  to  be  Im,  respectively,  and  the  standard 
deviation  a  is  5%.  The  backscattering  strength  versus  grazing  angle  can  be  estimated 
following  the  procedure  presented  in  the  beginning  of  this  section.  The  comparison 
with  the  analytic  solution  can  be  seen  in  Fig.  5-5,  in  which  the  error  bars  show  the 
standard  deviation  of  the  estimate.  To  our  delight,  the  numerical  solution  and  the 
analytic  solution  match  very  well. 

Second  is  an  example  of  backscattering  from  random  sound  speed  variations  with 
a  2-D  isotropic  exponential  correlation  function.  Figure  5-6  shows  the  ensemble  av¬ 
erage  over  50  realizations  of  simulated  backscattering  intensities  at  one  receiver.  All 
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Figure  5-5:  Comparison  of  the  backscattering  strength  estimated  from  the  numerically 
simulated  backscattered  returns  with  the  analytic  solution.  Sound  speed  variations 
have  a  Gaussian  correlation  function.  (/  =  400^/2,  H  =  400m,  =  L  =  Im,  a  = 

0.05) 


Figure  5-6:  Backscattering  from  2-D  isotropic  sound  speed  variations  with  an  ex¬ 
ponential  correlation  function,  averaging  over  50  realizations  (/  =  A00Hz,H  = 
400m,  Ij.  =  ly  =  0.84m,  a  =  0.05) 
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Figure  5-7:  Comparison  of  the  backscattering  strength  estimated  from  the  numeri¬ 
cally  simulated  backscattered  returns  with  the  analytic  solution.  Sound  speed  vari¬ 
ations  have  an  exponential  correlation  function.  (/  =  400Hz,H  =  400m,  L  =  L  = 
0.84m,  a  =  0.05) 


the  parameters  are  the  same  as  in  the  Gaussian  correlation  case  except  that  the  cor¬ 
relation  lengths  lx  and  ly  are  arbitrarily  chosen  to  be  0.84m,  respectively.  Again,  the 
backscattering  strength  versus  grazing  angle  is  estimated  and  the  comparison  with 
the  analytic  solution  is  plotted  in  Fig.  5-7.  The  matches  between  the  two  are  excellent. 

Third  is  an  example  of  backscattering  from  random  sound  speed  variations  with 
a  2-D  isotropic  power  law  distribution.  Figure  5-8  depicts  the  ensemble  average  over 
50  realizations  of  simulated  backscattering  intensities  at  one  receiver.  All  the  param¬ 
eters  are  the  same  as  in  the  Gaussian  correlation  case  and  the  horizontal  correlation 
lengths  are  once  again  Im  in  both  x  and  y  directions.  The  order  u  in  the  Von  Karman 
function  is  chosen  to  be  1  so  that  the  fractal  dimension  is  2.  This  is  the  case  for  the 
anisotropic  case  in  the  next  subsection  as  well.  Figure  5-9  shows  the  comparison  of 
the  estimated  backscattering  strength  with  the  analytic  solution.  The  results  again 
are  excellent. 
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Figure  5-8:  Backscattering  from  2-D  isotropic  sound  speed  variations  with  a  power 
law  distribution,  averaging  over  50  realizations  (/  =  400/fz,  H  =  400m,  lx  =  ly  = 
lm,a  =  0.05,  u  =  1) 


Figure  5-9:  Comparison  of  the  backscattering  strength  estimated  from  the  numerical 
simulated  backscattering  returns  with  the  analytic  solution.  Sound  speed  variations 
have  a  power  law  distribution.  (/  =  400Hz,  H  =  400m,  =  ly  =  Im,  a  =  0.05,  u  = 

1) 
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Figure  5-10:  Backscattering  from  2-D  anisotropic  sound  speed  variations  with  a 
Gaussian  distribution,  averaging  over  50  realizations  (  /  =  400Hz,  H  =  400m,  = 

Im,  ly  =  0.6m,  a  =  0.05) 

5.2.2  2-D  Anisotropic  Case 

The  numerical  model  will  be  tested  for  2-D  anisotropic  cases  in  this  subsection.  Again 
we  will  present  examples  for  three  different  types  of  correlation  functions. 

Figure  5-10  demonstrates  the  anisotropic  Gaussian  case.  The  plotted  backscat¬ 
tering  intensity  is  the  ensemble  average  of  50  realizations.  All  the  parameters  are 
the  same  as  in  the  isotropic  Gaussian  case  except  that  the  correlation  length  in  the 
y  direction  ly  has  been  changed  to  0.6m.  The  comparison  between  numerical  solu¬ 
tions  and  analytic  solutions  is  shown  in  Fig.  5-11.  For  an  anisotropic  exponential 
distribution,  the  backscattering  intensity  versus  arrival  time  is  depicted  in  Fig.  5-12. 
The  comparison  of  estimated  backscattering  strength  with  the  analytic  solution  can 
be  seen  in  Fig.  5-13.  Except  that  the  horizontal  correlation  length  ly  is  changed 
to  2m,  all  other  parameters  are  the  same  as  in  the  isotropic  Gaussian  case.  For  an 
anisotropic  power  law  distribution  with  correlation  length  l^  =  Im,  ly  =  0.6m,  and 
u  =  1,  Fig.  5-14  shows  the  backscattering  intensity  and  Fig.  5-15  the  comparison 
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Figure  5-11:  Comparison  of  the  backscattering  strength  estimated  from  the  numeri¬ 
cally  simulated  backscattered  returns  with  the  analytic  solution.  Sound  speed  varia¬ 
tions  have  a  Gaussian  distribution.  (/  =  AQQHz,  H  =  400m,  =  Im,  L  =  0.6m,  a  — 

0.05) 
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Figure  5-12.  Backscattering  from  2-D  anisotropic  sound  speed  variations  with  an 
exponential  distribution,  averaging  over  50  realizations  (/  =  400Hz,  H  =  400m,  /j,  = 
Im,  ly  =  2m,  a  =  0.05) 
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Figure  5-13:  Comparison  of  the  backscattering  strength  estimated  from  the  numerical 
simulated  backscattering  returns  with  the  analytic  solution.  Sound  speed  variations 
have  an  exponential  distribution.  (/  =  4QQHz,H  —  400m,  =  lm,ly  —  2m,  a  = 

0.05) 

with  analytic  solutions.  It  is  evident  that  matches  between  the  numerical  results 
and  the  analytic  solutions  are  excellent  for  all  three  of  the  above  anisotropic  cases. 


For  the  sake  of  brevity,  we  give  comparison  of  the  backscattering  strength  versus 
grazing  angle  only  for  one  set  of  parameters  in  each  distribution.  However,  we  have 
no  doubt  that  the  curves  in  Figs.  3-9,  3-10  and  3-11  can  be  reproduced  by  numerical 
simulation.  Since  the  decay  rate  of  the  backscattering  strength  versus  grazing  angle 
is  directly  related  to  the  decay  rate  of  backscattering  intensity  versus  time,  we  expect 
the  backscattering  intensity  to  have  different  decay  rates  versus  time  for  different 
values  of  the  parameters  in  the  above  three  figures.  Figure  5-16  gives  one  example  in 
the  Gaussian  distribution  case. 


Through  comparisons  with  analytic  solutions  in  2-D  isotropic  and  anisotropic  sit¬ 
uations,  the  monostatic  backscattering  model  has  been  proved  to  be  successful  in 
simulating  backscattered  returns  resulting  from  2-D  random  sound  speed  variations. 


Figure  5-14:  Backscattering  from  2-D  anisotropic  sound  speed  variations  with  a  power 
law  type  of  distribution,  averaging  over  50  realizations  (/  =  =  400m,  = 

Im,  ly  —  0.6m,  a  =  0.05,  v  =  1) 


Figure  5-15:  Comparison  of  the  backscattering  strength  estimated  from  the  numeri¬ 
cally  simulated  backscattered  returns  with  the  analytic  solution.  Sound  speed  varia¬ 
tions  have  a  power  law  distribution.  (/  =  AQQHz,  H  =  400m,  =  Im,  L  =  0.6m,  a  = 

0.05,  i/  =  1) 
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Figure  5-16:  Backscattering  intensities  for  simulated  backscattered  returns  due  to 
sound  speed  variations  with  a  Gaussian  distribution,  averaging  over  50  realizations 
(/  =  AmHz,  H  =  400m,  a  =  0.05). 


The  good  fit  between  the  numerical  results  and  analytic  solutions  also  implies  that 
the  1/8  wavelength  sampling  interval  is  sufficient  for  all  the  above  distributions. 


5.2.3  3-D  Vertically  Correlated  Case 

Three-dimensional  volume  scattering  with  vertical  correlation  is  more  complicated 
than  the  above  2-D  cases.  On  the  one  hand,  the  existence  of  physical  boundaries, 
i.e.,  the  layer  interfaces,  violates  the  stationarity  assumption  of  the  random  field, 
which  causes  the  edge  effect  that  does  not  arise  in  the  numerical  solution.  On  the 
other  hand,  arrival  time  can  no  longer  be  translated  into  scattering  angle  accurately. 
Returns  from  scatterers  located  at  different  depths  can  have  the  same  arrival  time, 
which  brings  in  angular  ambiguities  for  the  estimated  backscattering  strength  in  a 
single  source/single  receiver  configuration. 

In  this  section,  our  focus  is  on  thin  scattering  layers.  As  long  as  the  layer  thickness 
is  small,  the  single  source/single  receiver  configuration  is  still  capable  of  providing  a 
good  estimate  of  the  backscattering  strength  versus  grazing  angle.  However,  in  Chap¬ 
ter  2,  an  assumption  is  made  in  obtaining  the  analytic  solution  that  the  dimension  of 
the  scattering  volume  be  much  larger  than  the  correlation  length.  This  would  conflict 
with  the  small  layer  thickness  requirement  above.  A  compromise  has  to  be  reached 
in  choosing  the  layer  thickness. 

Again,  we  need  a  convergence  test  to  determine  the  sufficient  sampling  rate  in 
the  vertical  direction.  In  Fig.  5-17,  we  choose  a  7.5m  thick  random  layer  with  expo¬ 
nential  correlation  function  horizontally  and  vertically.  The  correlation  lengths  are 
lx  =  ly  =  2m  and  =  0.5m.  The  frequency  is  400Hz  and  the  source/receiver  height 
is  400m.  It  can  been  seen  that  the  curve  with  1/10  wavelength  as  sampling  interval 
in  the  2:  direction  is  very  close  to  that  with  1/20  wavelength.  Disregarding  the  be¬ 
ginning  part  which  corresponds  to  high  grazing  angles,  there  are  only  up  to  O.SdB 
discrepancies  along  the  curve  which  can  be  considered  acceptable.  Therefore,  1/10 
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Exponential  correlation  function  horizontally  and  vertically 


Figure  5-17:  The  convergence  test  for  backscattering  from  3-D  random  sound  speed 
variations  with  exponential  distributions  horizontally  and  vertically  (cr  =  0.05,  lx  = 
ly  —  2m,  Iz  =  0.5m,  /  =  400Hz,  H  =  400m,  D  =  7.5m) 

wavelength  is  regarded  as  a  sufficient  sampling  interval  in  the  vertical  direction. 

Following  the  same  procedure  as  in  the  2-D  cases,  the  backscattering  strength 
can  be  estimated  from  simulated  backscattered  returns.  The  comparison  with  the 
analytic  solution  is  shown  in  Fig.  5-18.  The  match  between  the  two  curves  is  not 
as  good  as  that  in  the  2-D  cases.  However,  except  for  regions  close  to  55  degree 
grazing  angle,  where  the  difference  is  about  2dB,  the  two  curves  are  close  to  each 
other.  Bear  in  mind  that  the  analytic  solution  is  obtained  under  the  assumption  that 
the  dimension  of  the  scattering  volume  is  much  larger  than  the  correlation  length  so 
that  the  boundary  effect  can  be  ignored,  while  the  numerical  solution  does  include 
the  boundary  effect.  That  could  well  be  the  reason  why  such  disagreement  exists. 
Anyway,  the  fit  can  be  regarded  as  satisfactory  since  the  numerical  solution  predicts 
the  trend  of  the  backscattering  strength  versus  grazing  angle  reasonably  well  and  the 
level  closely. 

In  conclusion,  the  numerical  solution  obtained  from  the  monostatic  backscattering 
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Figure  5-18;  Comparison  of  the  estimated  backscattering  strength  from  numerically 
simulated  backscattered  returns  with  the  analytic  solution.  Sound  speed  variations 
are  exponentially  correlated  horizontally  and  vertically,  {a  =  0.05,  =  ly  =  2m,  = 

0.5m,  /  =  400Hz,  H  =  400m,  D  =  7.5m) 

model  has  been  found  to  be  in  good  agreement  with  the  analytic  solution.  It  gives  us 
confidence  that  numerical  errors  are  under  control  and  the  model  prediction  is  valid, 
which  sets  the  stage  for  the  model/data  comparison  in  the  next  chapter. 


5.3  The  Effect  of  Density  Variations 

As  pointed  out  in  Chapter  3,  the  density  variations  act  like  a  dipole  source  with 
its  major  axis  in  the  backscattering  direction.  Therefore,  we  would  expect  total 
backscattering  to  be  enhanced.  To  confirm  this,  the  following  procedures  are  to  be 
taken:  First,  backscattering  returns  from  2-D  random  sound  speed  variations  only 
will  be  simulated.  Sound  speed  variations  are  assigned  a  power  law  distribution  with 
correlation  lengths  lx  —  ly  —  Im.  The  standard  deviation  is  5%.  The  frequency  is 
400Hz  and  the  source/receiver  height  is  400m.  The  backscattering  intensity  can  then 
be  obtained  and  the  backscattering  strength  estimated.  Next,  we  add  density  varia- 
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Figure  5-19:  Backscattering  intensities  for  simulated  backscattered  returns  due  to 
2-D  random  sound  speed  variations  only(dash-dotted  line)  and  density  variations 
only  (dash  line)  and  both(solid  line).  Sound  speed  variations  have  a  power  law  distri¬ 
bution.  (<T  =  0.05,  lx  =  ly  =  Im,  P  =  2,  f  =  AOOHz,  H  =  400m) 

tions  to  the  picture.  The  ratio  between  density  variations  and  sound  speed  variations, 
i.e.,  P  is  chosen  to  be  2  according  to  Hines’  empirical  results  [14].  The  backscattering 
intensity  is  recalculated  and  the  backscattering  strength  again  estimated. 

Backscattering  intensities  with  sound  speed  variations  only,  with  density  varia¬ 
tions  only  and  with  both  variations  are  plotted  together  in  Fig.  5-19  and  the  cor¬ 
responding  backscattering  strengths  in  Fig.  5-20.  The  backscattering  strength  with 
density  variations  only  is  7dB  higher  than  that  with  sound  speed  variations  only, 
which  shows  that  density  variations  are  the  more  dominant  contributors  to  backscat¬ 
tering.  An  interesting  phenomenon  is  that  the  two  curves  are  parallel  to  each  other.  In 
other  words,  the  angular  dependence  of  the  backscattering  strength  has  not  changed 
with  the  addition  of  density  variations.  One  possible  explanation  is  that  the  major 
axis  of  the  virtual  dipole-like  source  caused  by  density  variations  always  lies  in  the 
backscattering  direction.  While  density  variations  and  sound  speed  variations  are 
fully  correlated  in  the  meantime,  the  dipole-like  source  will  behave  the  same  as  the 
monopole-like  source  caused  by  sound  speed  variations  in  the  backscattered  direction. 
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Figure  5-20:  Backscattering  strengths  for  simulated  backscattered  returns  due  to 
2-D  random  sound  speed  variations  only(dash-dotted  line)  and  density  variations 
only  (dash  line)  and  both  (solid  line)  .  Sound  speed  variations  have  a  power  law  distri¬ 
bution.  (a  =  0.05,  lx  =  ly  =  im,l3  =  2,f  =  400Hz,  H  =  400m) 


All  in  all,  the  above  examples  have  clearly  demonstrated  the  importance  of  in¬ 
cluding  density  variations  in  backscattering  modeling.  Since  density  variations  will 
change  only  the  level  but  not  the  angular  dependence  of  the  backscattering  strength, 
the  conclusions  that  we  have  reached  through  the  study  of  sound  speed  variations 
only  will  still  be  valid. 
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Chapter  6 


Comparison  of  Model  Results  and 
ARSRP  Backscattering  Data 

A  primary  goal  of  this  work  is  to  develop  a  volume  backscattering  model  that  can 
account  for  the  backscattering  phenomena  observed  in  the  ARSRP  sediment  experi¬ 
ment.  Throughout  Chapter  2,  3,  and  4,  the  accumulating  efforts  have  clearly  paved 
the  way  for  the  realization  of  a  backscattering  model,  which  is  capable  of  incorporat¬ 
ing  contributions  due  to  both  sound  speed  and  density  variations.  In  the  meanwhile, 
by  comparing  to  analytic  solutions  in  free-space  scenarios,  the  reliability  of  the  model 
has  been  demonstrated  in  the  previous  chapter.  Nevertheless,  a  more  important  but 
challenging  task  is  to  evaluate  the  capability  of  the  model  in  predicting  the  experi¬ 
mental  results,  which  is  the  concentration  of  this  chapter. 

As  discussed  earlier,  an  ideal  and  straightforward  way  to  compare  the  model  and 
backscattering  data  due  to  volume  scattering  from  a  thick  layer  of  sediment  inhomo¬ 
geneities  is  through  the  scattering  time  series.  However,  in  the  ARSRP  experiment 
scenario,  the  wide  bandwidth  of  the  source  signal,  400R^,  together  with  the  large 
layer  thickness,  20  —  30m  (has  to  be  sampled  at  every  1/lOth  wavelength  vertically 
in  order  to  guarantee  convergence)  make  the  simulation  of  the  scattering  time  series 
extremely  computationally  intensive.  It  would  be  unrealistic  in  this  situation  to  run 
the  model  many  times  so  as  to  achieve  a  good  fit  with  the  data.  Therefore,  an  al- 
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ternative  scheme  for  comparison  is  chosen,  which  is  to  compare  the  backscattering 
strength  versus  grazing  angle  between  model  and  data  for  selected  frequencies  across 
the  band.  One  might  recall  that  the  single  source/single  receiver  configuration  cannot 
distinguish  scattered  returns  from  different  directions  but  with  the  same  arrival  time, 
thereby  resulting  in  angular  ambiguities  of  the  backscattering  strength.  Fortunately 
for  us,  this  drawback  can  be  overcome  thanks  to  the  vertical  line  array  and  beam¬ 
forming  technique.  Besides,  as  long  as  the  model  and  data  are  processed  by  the  same 
means,  the  results  of  comparisons  will  not  be  affected. 

This  chapter  begins  by  reviewing  the  experimental  setup.  All  the  information 
about  the  backscattering  data  collection  is  provided  in  the  first  section.  The  ARSRP 
data  processing  will  then  be  briefly  discussed  in  the  second  section.  In  the  third 
section,  the  focus  will  be  on  the  estimation  of  the  backscattering  strength  from  the 
beamformed  backscattered  returns.  Details  will  also  be  given  on  how  to  choose  the 
model  parameters  since  they  are  closely  related  to  the  above  estimation.  The  com¬ 
parison  between  model  results  and  the  ARSRP  backscattering  data  is  the  theme  of 
the  final  section. 


6.1  Experiment  Description 

On  July  24-26, 1993,  as  part  of  the  Acoustic  Reverberation  Special  Research  Program 
(ARSRP),  a  bottom  backscattering  experiment  was  conducted  over  the  site  A  sedi¬ 
ment  pond  located  at  26°11’N,  46°09’W  in  the  vicinity  of  the  Mid-Atlantic  Ridge  [76]. 

The  water  depth  was  about  4400  meters  for  most  of  the  sediment  pond  area  and 
the  sound  speed  profile  was  a  typical  linear  deep  ocean  profile  obtained  from  XBT 
measurements.  The  gradient  was  0.0154/s  for  the  water  column  close  to  the  seafloor. 
The  sound  speed  at  the  water-sediment  interface  was  extrapolated  to  be  1530  m/s. 
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The  oval-shaped  sediment  pond  was  approximately  15  km  x  20  km  in  size.  The 
thickness  of  the  sediment  cover  varied  from  almost  no  cover  at  the  edge  of  the  pond 
to  about  430  meters  in  the  middle  of  the  pond.  For  areas  where  the  backscattered 
data  were  collected,  the  sediment  was  at  least  100  m  thick.  The  top  layer  beneath  the 
water-sediment  interface  was  made  of  calcareous  mud  with  a  sound  speed  of  1510  m/s 
according  to  core  data  [77].  No  direct  measurements  of  the  attenuation  coefficient  or 
the  sound  speed  gradient  were  available  in  the  sediment  area  where  the  experiment 
was  conducted.  A  deep  towed  high-frequency  chirp  sonar  and  in  situ  optical  survey 
revealed  that  the  seafloor  was  virtually  flat  with  rms  roughness  less  than  10  cm  [78]. 
A  rough  basalt  basement  was  beneath  the  sediment  layer. 

Figure  6-1  shows  the  bathymetry  of  the  sediment  pond  and  the  ship  track  along 
which  the  experiment  was  carried  out.  For  the  results  presented  in  this  chapter,  we 
will  only  analyze  the  data  collected  along  the  straight  line  shown  in  the  figure.  Fig¬ 
ure  6-2  depicts  the  experimental  scenario.  The  acoustic  system  was  deployed  at  the 
west  side  of  the  pond  and  recovered  in  the  middle  after  being  towed  across  the  entire 
sediment  pond. 

The  acoustic  source  used  in  the  experiment  is  part  of  the  Deep  Towed  Acous¬ 
tics/Geophysics  System  (DTAGS)  and  was  developed  by  the  Naval  Research  Labo¬ 
ratory  -  Stennis  Space  Center [79].  The  source  signal  is  a  linear  frequency-modulated 
(LFM)  upchirp  signal  in  the  frequency  range  250-650  Hz.  The  designed  source  level  is 
200  dB  re  1//Pa  @  1  m.  The  duration  of  the  source  signal  is  0.125  s.  The  receiver  is  a 
24  element  vertical  line  array(VLA)  hung  beneath  the  source,  with  a  spacing  between 
adjacent  hydrophone  elements  of  2  m.  The  cable  connecting  the  source  and  the  top 
hydrophone  is  57.5  m  long.  A  weight  is  attached  to  the  end  of  the  array  to  keep 
it  close  to  a  vertical  configuration.  Figure  6-3  shows  the  geometry  of  the  acoustic 
transmitting  and  receiving  system.  The  source  signal  is  transmitted  once  every  32 
s.  For  each  ping,  the  hydrophones  record  4.5  s  of  waveforms  starting  0.1  s  before 
the  source  trigger.  The  sampling  interval  is  0.432  ms.  More  than  5000  pings  were 
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Figure  6-1:  Bathymetry  and  ship  track  in  the  ARSRP  Site  A  sediment  pond. 
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Figure  6-2:  The  experimental  scenario. 
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Figure  6-3:  The  DTAGS  source  and  receiving  array  geometry. 

transmitted  throughout  the  experiment.  Figure  6-4  shows  the  received  signal  and  its 
spectrum  at  one  hydrophone  for  one  ping. 

During  the  experiment,  the  DTAGS  source  was  towed  at  a  depth  of  about  200- 
400  meters  above  the  seafloor,  while  the  ship  used  the  high-resolution  (P-code)  GPS 
dynamic  positioning  to  maintain  a  constant  speed  over  ground  (SOG)  of  0.5  knot, 
with  several  stops  to  ensure  that  the  array  returned  to  nearly  vertical  alignment.  The 
actual  array  shape  is  estimated  by  using  the  direct  and  bottom  reflected  arrival  times 
on  each  hydrophone.  The  details  of  the  array  element  localization  method  was  given 
in  [10].  Figure  6-5  shows  the  estimated  array  configuration  for  one  location.  It  can 
be  noticed  that  the  array  is  roughly  vertical.  Generally  it  is  found  that  the  array  tilt 
angle  is  within  5  degrees  of  vertical. 
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Figure  6-4:  (a)  The  received  signal  and  (b)  its  spectrum  at  one  hydrophone  for  a 
typical  ping 
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Figure  6-5:  The  estimated  receiving  array  configuration  for  a  typical  ping 
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6.2  Data  Processing 


Only  a  short  version  of  the  data  processing  techniques  will  be  presented  here  since 
details  have  been  given  in  Ref.  [10].  The  emphasis  of  this  section  is  on  estimating 
the  backscattering  strength  versus  grazing  angle  from  the  data. 

The  sediment  is  well-known  for  its  layered  structure,  which  is  further  confirmed 
by  bottom  profiling  shown  in  Figs.  1-1  and  1-3  in  the  experimental  region.  It  is  evi¬ 
dent  from  the  plot  that  strong  coherent  reflections  exist  at  layer  interfaces.  With  the 
sediment  thickness  ranging  from  100m  to  400m,  the  reflections  from  deep  layers  can 
arrive  at  the  same  time  as  backscattered  returns  in  oblique  directions.  As  a  result, 
the  usual  time-gating  technique  are  no  longer  valid  in  separating  the  reflected  and 
backscattered  returns.  In  other  words,  the  assumption  that  scattering  arrivals  at  a 
certain  time  are  associated  with  only  one  scattering  angle  are  not  valid.  This  can 
explain  why  the  angular  dependence  of  the  backscattering  strength  obtained  from  a 
single  source/single  receiver  configuration  is  error-prone. 

The  aforementioned  problem  had  been  recognized  in  the  design  of  the  ARSRP 
backscattering  experiment  and  a  vertical  line  array  instead  of  a  single  hydrophone 
was  used  as  the  receiver.  With  the  beamforming  technique,  scattered  returns  in  the 
look  direction  can  be  obtained  without  distortion,  while  all  the  returns  in  other  di¬ 
rections  have  been  attenuated.  The  beam  pattern  for  a  conventional  single-constraint 
beamformer  is  shown  in  Fig.  6-6.  The  array  here  has  the  same  geometry  as  the  actual 
receiving  array  depicted  in  Fig.  6-3,  i.e.,  a  24  element  vertical  line  array  with  2m 
spacing  between  hydrophones.  The  frequency  is  chosen  to  be  the  center  frequency  of 
the  source  signal,  450/^2;  and  the  look  direction  is  — 45  degrees.  The  array  response 
for  each  circle  is  20dB  higher  than  that  of  the  neighboring  inner  circle.  The  grating 
lobes  in  the  upper  part  of  the  plot  are  caused  by  undersampling.  From  the  plot¬ 
ted  beam  pattern,  as  expected,  one  can  see  that  no  attenuation  is  imposed  on  signals 
coming  from  the  look  direction.  However,  in  sediment  volume  scattering,  the  coherent 
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Figure  6-6:  The  beam  pattern  for  a  single-constraint  beamformer:  frequency  =  450 
Hz,  look  direction  =  —45  degrees. 
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Figure  6-7:  The  beam  pattern  for  a  five-constraints  beamformer:  frequency  =  450 
Hz,  look  direction  =  —45  degrees 

reflections  in  the  normal  and  near  normal  incidence  directions  are  much  stronger  than 
the  backscattered  returns  in  oblique  directions,  which  had  been  observed  in  ARSRP 
backscattering  experiment.  Therefore,  there  was  a  great  possibility  that  the  above 
shown  sidelobe  levels  in  the  normal  and  near  normal  incidence  directions  were  not 
low  enough  to  reject  those  unwanted  returns.  The  ineffectiveness  of  the  conventional 
beamformer  had  been  confirmed  by  simulations  [10].  To  overcome  this,  a  multiple- 
constraints  beamformer  was  developed  [10]  and  its  beampattern  for  the  same  setup  as 
in  the  above  conventional  single-constraint  beamformer  is  shown  in  Fig.  6-7.  A  null 
is  formed  in  the  normal  and  near  normal  incidence  directions,  thereby  improve  the 
beamformer’s  ability  of  alleviating  the  strong  reflections.  Simulations  in  Ref.  [10]  saw 
much  more  favorable  results  compared  to  those  with  conventional  single-constraint 
beamforming. 
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To  further  reduce  the  normal-incidence  returns,  a  hybrid  method,  which  com¬ 
bines  multiple  constraints  beamforming  and  a  subtraction  procedure,  was  applied  in 
the  ARSRP  data  processing.  The  idea  is  to  subtract  first  the  strong  returns  from 
the  normal  and  near  normal  incidence  directions  obtained  from  endfire  beamforming, 
then  apply  multiple  constraints  beamforming  to  estimate  the  backscattered  returns 
in  oblique  directions.  Bear  in  mind  that  it  is  case-dependent  for  the  subtraction  pro¬ 
cedure  to  work  and  it  is  almost  impossible  to  subtract  all  the  returns  in  the  normal 
and  near  normal  incidence  directions  for  experimental  data.  However,  the  subtrac¬ 
tion  proved  to  be  successful  for  the  ARSRP  backscattering  data,  since  most  of  the 
high  amplitude  peaks  due  to  reflections  disappeared  after  the  subtraction  (cf.  Fig. 
3-11  [10]).  The  multiple  constraints  beamformer  will  help  to  reject  the  residue  of  the 
subtraction,  which  is  also  in  the  normal  and  near  normal  incidence  directions,  and 
obtain  the  uncontaminated  backscattered  returns  in  oblique  directions.  This  led  to 
the  findings  stated  in  Chapter  1  that  two  irregular  layers  beneath  the  water/sediment 
interface  are  the  primary  contributors  to  backscattering  in  oblique  directions. 


6.3  Selection  of  Model  Parameters  and  Estimation 
of  the  Backscattering  Strength 

Since  the  scattering  layers  are  beneath  the  seafloor,  the  selection  of  background  pa¬ 
rameters  for  the  model  such  as  the  sound  speed  profile  and  attenuation  coefficient 
would  directly  affect  the  estimation  of  the  backscattering  strength,  which  prompts  us 
to  address  these  two  together  in  this  section. 

The  extrapolation  of  the  CTD  measurement  revealed  that  the  water  sound  speed 
at  the  top  of  the  water/sediment  interface  is  1530m/s.  According  to  Fu  et.  al.  [77], 
the  sound  speed  just  beneath  the  interface  is  1510m/s,  slightly  lower  than  that  in 
the  water  column.  Using  a  ray-tracing  type  of  estimation  of  the  sediment/basement 
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Figure  6-8:  Fluid  bottom  reverberation  scenario. 


interface  returns,  the  average  sound  speed  gradient  in  the  sediment  is  found  to  be 
approximately  1.12/s  [9].  According  to  Hamilton  [62],  a  typical  gradient  for  this  type 
of  sediment  is  about  1/s.  From  the  amplitudes  of  the  normal-incidence  reflections  at 
the  water /sediment  interface  and  the  sediment/basement  interface  together  with  the 
arrival  time,  the  attenuation  is  approximately  estimated  as  0.\dB/\,  which  is  consis¬ 
tent  with  values  in  the  literature  [59,  32]  for  the  calcareous  mud  type  of  sediment  in 
this  area.  The  density  of  the  sediment  is  chosen  to  be  ISOOA:^/?^^,  also  well  within 
the  range  for  this  type  of  sediment.  By  incorporating  all  the  above  parameters,  we 
propose  a  simple  fluid  bottom  model  depicted  in  Fig.  6-8.  The  sound  speed  is  as¬ 
sumed  to  be  linear  from  the  water /sediment  interface  down  to  100m.  The  interface  at 
100m  beneath  the  water /sediment  interface  is  a  false  one,  as  is  the  isovelocity  bottom 
beneath  it.  As  a  result,  the  specular  reflections  at  the  sediment/basement  interface 
and  their  possible  insonification  of  the  two  irregular  scattering  layers  are  neglected 
because  they  would  arrive  later  than  backscattered  returns  at  30  degree  grazing  angle. 
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which  is  the  smallest  angle  at  which  backscattered  returns  can  still  be  distinguished 
from  the  background  noise  (cf.  Fig.  1-2).  The  positions  of  the  two  irregular  layers 
can  be  determined  from  the  travel  time  obtained  from  bottom  profiling  shown  in  Fig. 
1-1  together  with  the  above  proposed  sound  speed  profile.  A  simple  calculation  indi¬ 
cates  that  the  the  first  irregular  layer  is  located  from  H\  =  15.85m  to  H2  —  31.42m 
beneath  the  water /sediment  interface,  while  the  second  one  is  from  Hz  =  56.42m 
to  Hi  =  85.7m.  Notice  that  these  numbers  are  different  from  those  marked  on  the 
side  of  the  bottom  profiling  plot  because  a  nominal  sound  speed  of  1530m/s  was  used 
then.  The  water/sediment  interface  was  found  to  be  flat  as  mentioned  in  the  previous 
section  and  verified  by  the  data  analysis  in  [10,  13],  which  vindicates  the  choice  of  a 
range-independent  water/sediment  interface  with  no  roughness  in  the  model.  All  the 
above  parameters  will  be  the  input  to  OASES  so  that  the  Green’s  function  can  be 
evaluated. 

To  characterize  the  backscattering  process  quantitatively,  we  still  resort  to  estima¬ 
tions  of  the  backscattering  strength.  As  you  can  recall  in  Chapter  2,  we  had  discussed 
the  need  to  be  cautious  in  using  the  backscattering  strength  to  describe  a  volume 
scattering  process.  This  quantity  might  be  inappropriate  when  multi-path  effects  are 
significant.  However,  the  following  reasons  make  us  believe  that  the  backscattering 
strength  versus  grazing  angle  is  an  acceptable  parameter  in  describing  backscattering 
from  the  two  irregular  layers  found  in  the  ARSRP  experiment:  First,  thanks  to  the 
vertical  line  array  and  beamforming,  we  can  discriminate  the  oblique  backscattered 
returns  from  reflections  at  layer  interfaces,  which  makes  the  angular  dependence  of 
the  estimated  backscattering  strength  trustworthy;  second,  the  multi-path  effect  is 
not  significant  because  of  a  lower  sound  speed  at  the  top  of  the  sediment  compared 
to  that  in  the  water  right  above  it  and  because  of  the  small  sound  speed  gradient 
in  the  bottom;  in  addition,  we  are  only  interested  in  the  angular  regime  between  30 
to  60  degrees  grazing  angles,  where  the  sediment/basement  reflections  arrive  late  as 
discussed  before. 
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The  estimation  of  the  backscattering  strength  from  data  is  no  simple  task.  Al¬ 
though  the  procedures  were  described  in  Ref.  [10]  and  reiterated  in  Chapter  5,  we 
would  like  to  give  more  details,  particularly  on  how  to  arrive  at  the  frequency  depen¬ 
dence  of  the  backscattering  strength.  Let’s  start  from  Eq.5.1 

BSS  =  10  log,„  0, 

where  is  the  intensity  of  the  backscattered  signals  at  the  reference  receiver,  /i„  is  the 
intensity  of  the  incident  sound  wave  a  unit  distance  away  from  the  scattering  volume, 
R  is  the  distance  from  the  reference  receiver  to  the  top  of  the  scattering  volume,  and 
A  is  the  top  area  of  the  scattering  volume.  A  reference  receiver  is  needed  because  we 
now  have  a  vertical  line  array  instead  of  a  single  receiver.  In  this  data  analysis,  the 
middle  hydrophone  of  the  array  (the  twelfth  from  the  top)  is  chosen  as  the  reference 
receiver.  Also,  the  directivity  aspect  of  the  array  and  subsequent  beamforming  makes 
it  difficult  to  determine  the  insonified  area  precisely.  Here,  we  just  consider  the  simple 
case  where  the  returns  are  from  the  maximum  response  axis(MRA)  of  the  array,  i.e., 
the  returned  signals  suffer  little  attenuation  passing  through  the  beamformer.  As  a 
result,  A  is  estimated  by  the  3dB  area  corresponding  to  the  annulus  where  the  array 
response  is  within  SdB  compared  to  the  MRA.  By  applying  ray  tracing,  A  can  be 
calculated  by 


A  =  7r(r|  -  rl),  (6.1) 

where  ri  and  r2  are  the  horizontal  ranges  between  the  reference  receiver  and  the 
radial  bounds  of  the  SdB  area(annulus)  at  the  top  of  the  scattering  layer,  as  shown 
in  Fig.  6-9.  Bear  in  mind  that  the  SdB  beamwidth  Oscib  and  therefore  the  SdB 
insonified  area  will  change  with  frequency.  The  lower  the  frequency,  the  wider  the 
beam,  the  larger  the  insonified  area.  In  order  to  get  the  backscattering  strength 
for  each  frequency,  we  need  to  estimate  Iin{f)  and  /«(/),  where  /  stands  for  the 
frequency.  To  calculate  /{„(/),  the  source  signal  is  first  Fourier  transformed.  The 
amplitude  squared  of  the  source  spectrum  is  then  multiplied  by  R?  to  compensate 
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for  the  spreading  loss.  Meanwhile,  the  calculation  of  Is{f)  would  involve  one  more 
procedure,  which  is  to  determine  the  exact  piece  of  beamformed  data  at  each  angle 
that  corresponds  to  the  backscattered  returns  from  each  scattering  layer.  It  can  be 
expressed  as  follows, 

/,(/)=<  IffTlyC*)}!"  >,  T,  <t<T2  (6.2) 

where  y{t)  are  the  scattered  returns,  FFT{-}  stands  for  the  Fourier  transform,  |  •  | 
stands  for  modulus,  <  •  >  stands  for  ensemble  average  over  different  pings,  and  Ti 
and  T2  are  the  travel  times  corresponding  to  arrivals  from  the  top  and  the  bottom  of  a 
scattering  layer.  Note  that  y{t)  along  with  Ti  and  T2  are  different  for  each  scattering 
angle.  The  backscattering  strength  versus  grazing  angle  for  each  frequency  can  then 
be  evaluated  according  to  Eq.5.1.  The  two  irregular  scattering  layers  are  treated  sep¬ 
arately.  When  computing  the  incident  field  on  the  lower  irregular  layer,  the  effect  of 
the  upper  irregular  layer  on  the  transmitted  field  beneath  it  is  neglected.  Actually  it  is 
consistent  with  the  Born  approximation  employed  in  our  scattering  model  since  scat¬ 
tering  from  volume  inhomogeneities  is  supposed  to  be  weak.  In  Tang’s  thesis  [35], 
simulations  had  demonstrated  that  the  transmitted  field  beneath  a  six-wavelength 
thick  random  layer  calculated  by  the  Born  approximation  was  accurate  enough.  Co¬ 
incidentally,  the  upper  irregular  layer  in  the  ARSRP  experiment  has  a  thickness  of 
about  six  wavelengths.  Therefore,  the  neglect  of  the  influence  on  the  incident  field 
of  the  lower  irregular  layer  due  to  the  existence  of  the  upper  irregular  layer  can  be 
justified.  Taking  into  consideration  the  attenuation  in  the  sediment,  an  amount  of 
2aRi  will  be  added  to  the  above  backscattering  strength,  where  a  is  the  attenuation 
coefficient  and  Ri  is  the  ray-path  length  in  the  sediment. 

Figures  6-10  and  6-11  show  the  estimated  backscattering  strength  versus  grazing 
angle  and  frequency  for  the  two  irregular  scattering  layers  at  the  east  side  of  the 
sediment  pond.  Note  that  the  grazing  angle  here  and  from  now  is  the  angle  measured 
at  the  water /sediment  interface  instead  of  at  the  top  of  the  irregular  region.  How- 
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Figure  6-10:  The  estimated  backscattering  strength  for  the  upper  irregular  region  at 
the  east  side  of  the  sediment  pond. 


Estimated  from  ARSRP  Data 


Figure  6-11:  The  estimated  backscattering  strength  for  the  lower  irregular  region  at 
the  east  side  of  the  sediment  pond. 
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ever,  due  to  the  fact  that  the  sound  speed  gradient  is  small,  the  largest  difference, 
appearing  at  30  degrees,  is  1.4  degrees  for  the  upper  layer  and  3.4  degrees  for  the 
lower  layer,  values  which  are  smaller  than  the  beamwidth.  Therefore,  the  angular 
dependence  obtained  can  be  considered  as  reliable.  It  can  be  seen  that  in  general 
the  backscattering  strength  would  decay  with  the  increase  of  the  frequency  and/or 
the  decrease  of  the  scattering  grazing  angle.  The  widely  used  Lambert’s  law  will  cer¬ 
tainly  fail  to  produce  the  frequency  dependence  of  the  backscattering  strength.  With 
a  point  scatter  model,  the  predicted  backscattering  strength  would  be  independent 
of  scattering  angle,  which  is  also  not  the  case  here.  In  Fig.  6-11,  there  is  a  peak 
standing  out  at  low  frequencies  and  high  grazing  angles.  The  peak  is  only  confined  to 
grazing  angles  larger  than  50  degrees  and  frequencies  less  than  350i/^,  which  makes 
it  look  more  like  an  artifact.  One  possible  cause  of  it  could  be  as  follows:  Although 
the  aforementioned  subtraction  procedure  had  eliminated  some  strong  normal  and 
near  normal  incidence  reflections,  it  can  hardly  be  perfect  in  a  real  environment.  As 
a  result,  there  will  be  some  residue  or  even  some  error  brought  in.  Certainly  the 
tilting  of  the  layers  beneath  the  lower  irregular  layer,  as  can  be  seen  in  Fig.  1-1, 
makes  matters  worse.  The  subtraction  process  would  face  greater  difficulties  here. 
Meanwhile,  the  beamwidth  of  the  mainlobe  is  larger  for  a  lower  frequency  and  when 
the  look  direction  is  closer  to  normal  incidence,  which  increases  the  chance  of  picking 
up  those  residues  in  the  normal  and  near  normal  incidence  directions.  Incidentally, 
the  returns  from  the  above  tilting  layers  would  have  about  the  same  arrival  time  as 
the  scattered  returns  from  the  lower  irregular  layers  at  about  a  60  degree  grazing 
angle.  In  addition,  ambient  noise  is  found  to  peak  between  250Hz  and  SOOHz,  which 
may  influence  the  scattering  strength  estimation  because  of  weak  oblique  backscat¬ 
tering.  Therefore,  in  the  model/data  comparisons,  we  will  concentrate  more  on  the 
backscattering  strength  at  frequencies  larger  than  300Hz. 

Notice  that  the  estimation  of  the  backscattering  strength  from  data  is  related  to 
the  background  bottom  model,  especially  the  attenuation  coefficient.  Because  the 
scattered  returns  from  the  two  irregular  layers  are  identified  by  their  travel  times. 
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Figure  6-12:  The  sensitivity  of  the  estimated  backscattering  strength  to  the  attenua¬ 
tion  coefficient  at  A^QHz  for  the  upper  irregular  layer  at  the  east  side  of  the  sediment 
pond. 

errors  in  the  sound  speed  profile  would  have  little  impact  on  the  backscattering  es¬ 
timation.  However,  the  attenuation  coefficient  would  decide  the  compensation  level 
(cf.  pl55)  of  the  backscattering  strength  at  different  scattering  angles,  therefore  in¬ 
fluencing  the  angular  dependence  of  the  backscattering  strength  as  well  as  the  level. 
In  this  study,  we  obtain  the  best  estimate  of  the  attenuation  coefficient  from  the 
normal  incidence  returns  as  According  to  the  archival  data  compiled  by 

Kibblewhite  [59],  the  attenuation  coefficient  for  silt  and  clay  sediment  is  between 
0.01d5/A  to  0.12dS/A.  The  example  shown  in  Fig.  6-12  illustrates  the  sensitivity  of 
the  estimated  backscattering  strength  at  450i^^  to  the  attenuation  coefficient. 
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6.4  Model  and  Data  Comparison 


By  no  means  it  is  an  easy  task  computationally  to  generate  scattering  time  series,  es¬ 
pecially  the  wide  400i72:  bandwidth  in  the  ARSRP  experiment.  Adding  to  the  degree 
of  difficulty  are  the  thick  scattering  layer  and  fine  sampling  requirement.  In  order 
to  perform  beamforming,  scattered  returns  have  to  be  simulated  across  the  whole 
receiving  array.  Therefore,  instead  of  a  scattering  time  series,  we  choose  to  compare 
the  backscattering  strength  versus  grazing  angle  estimated  from  model  simulation 
and  real  data  for  a  selected  set  of  frequencies  over  the  signal  spectrum. 

The  fact  that  multiple  constraints  beamforming  is  carried  out  in  the  frequency 
domain  provides  us  with  extra  benefit.  In  other  words,  if  the  scattered  field  at  each 
hydrophone  for  a  single  frequency  is  available,  we  are  still  able  to  perform  beam¬ 
forming  to  obtain  the  scattered  field  and  therefore  the  scattered  intensity  at  different 
scattering  angles.  As  a  result,  in  order  to  compare  with  data,  all  we  need  to  do  in 
the  modeling  part  is  to  simulate  the  scattered  field  at  each  receiver  for  the  selected 
set  of  frequencies.  The  beamforming  is  the  same  as  that  in  the  data  processing.  So 
is  the  estimation  of  the  backscattering  strength  after  obtaining  the  backscattering 
intensity  for  each  frequency.  Meanwhile,  the  procedures  of  simulating  the  scattering 
time  series  have  been  given  in  Chapter  5.  We  simply  forego  the  last  step,  the  Fourier 
synthesis,  since  our  goal  here  is  to  obtain  the  scattered  field  at  certain  frequencies. 

Model/data  comparisons  will  be  carried  out  for  the  two  irregular  layers  at  the  east 
side  of  the  sediment  pond  only.  We  choose  to  compare  at  six  frequencies:  300/fz, 
SbOHz,  400Hz,  4bOHz,  500Hz,  and  550Hz.  As  for  the  power  spectrum  and  char¬ 
acteristic  length  scale  to  describe  the  inhomogeneous  sound  speed  and  density  fluc¬ 
tuations,  we  elect  to  test  the  Gaussian  distribution  and  the  power  law  distribution 
(the  exponential  distribution  is  a  special  case  of  a  power  law  distribution).  There  is 
no  evidence  suggesting  horizontal  anisotropy  in  the  ARSRP  Site  A  sediment  pond 
area.  So  we  will  simply  let  the  horizontal  correlation  length  =  ly  in  the  model. 
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Meanwhile,  there  is  no  measurement  of  the  correlation  length  in  this  sediment  area. 
By  examining  Fig.  1-1,  one  would  have  the  impression  that  the  random  field  is  much 
better  correlated  horizontally  than  vertically.  As  a  matter  of  fact,  this  is  consistent 
with  the  geoacoustical  model  proposed  by  Lysanov  [29],  in  which  he  suggested  large- 
scale  in  the  horizontal  plane  and  small-scale  in  depth  for  random  inhomogeneities  in 
the  sediment,  the  so-called  “pancake”  model.  The  above  sediment  profiling  covered 
a  horizontal  distance  of  about  330m.  The  rough  estimation  would  show  us  that  the 
horizontal  correlation  length  be  on  the  order  of  10m.  Of  course,  the  estimation  can 
be  very  wrong.  However,  it  is  just  the  starting  point  of  our  parameter  search.  In  the 
meantime,  the  standard  deviation  o  of  the  sound  speed  fluctuation  can  range  from 
1.5%  to  8%  and  /?,  the  ratio  between  density  and  sound  speed  fluctuations,  between 
1  to  10  according  to  Yamamoto  [64].  Also  in  Yamamoto’s  estimation,  both  sound 
speed  and  density  fluctuations  were  characterized  by  a  power  law  distribution  with  v 
(cf.  Eq.3.35)  between  0.4  and  0.65. 

As  shown  in  Fig.  6-13,  the  best  fit  for  the  upper  irregular  layer  is  with  a  power  law 
distribution.  The  parameters  for  the  model  are  correlation  lengths  =  ly  =  20m, 
Zz  =  0.7m,  the  standard  deviation  of  sound  speed  fluctuations  a  =  2.8%,  j3  =  3, 
and  1/  =  0.5.  The  model  predicted  curves  are  ensemble  averages  of  200  realizations, 
therefore  the  standard  deviations  are  very  small,  (see  Appendix  B  for  the  estima¬ 
tion  of  the  standard  deviation  from  the  simulation  and  data.)  In  general,  the  model 
curves  agree  with  the  data  very  well.  Notice  that  in  the  data  curves,  around  a  47 
degree  grazing  angle  for  all  the  selected  frequencies  except  500Hz,  there  is  a  small 
bump.  No  specific  cause  has  been  found  for  this  and  the  model  cannot  predict  it 
either.  As  for  the  parallel  shift  between  the  model  and  data  curves  at  300Hz,  it 
might  be  due  to  the  aforementioned  strong  ambient  noise.  Figure  6-14  shows  the 
model/data  comparison  for  the  lower  irregular  layer  at  the  east  side  of  the  sediment 
pond.  The  best  fit  again  belongs  to  random  sound  speed  fluctuations  with  a  power 
law  type  of  power  spectrum.  The  horizontal  correlation  lengths  are  =  ly  =  20m, 
the  vertical  correlation  length  is  0.8m,  the  standard  deviation  a  =  2%,  P  =  3,  and 


161 


Backscattering  strength(dB) 


f=300Hz 


f=350Hz 


fi^OOHz  50 


60 


fi^SOHz  50  60 


Figure  6-13:  Model/data  comparison  for  the  upper  irregular  layer  at  the  east  side  of 
the  sediment  pond.  Sound  speed  fluctuations  are  described  by  a  power  law  distribu¬ 
tion  with  Ij.  —  ly  =  20m,  =  0.7m,  a  =  2.8%,  /?  =  3,  and  u  =  0.5.  The  error  bars 

show  the  standard  deviation  of  the  simulated  backscattering  strength. 
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Figure  6-14:  Model/data  comparison  for  the  lower  irregular  layer  at  the  east  side  of  the 
sediment  pond.  Sound  speed  fluctuations  are  described  by  a  power  law  distribution 
with  =  ly  =  20m,  =  0.8m,  cr  =  2%,  =  3,  and  i/  =  0.5.  The  error  bars  show  the 

standard  deviation  of  the  simulated  backscattering  strength. 


ly  =  0.5.  200  realizations  are  used  to  obtain  the  ensemble  average.  The  matches  are 
again  satisfactory,  although  the  error  in  sound  speed  profile  and  attenuation  estima¬ 
tion  would  have  more  impact  here  than  on  the  comparison  of  the  upper  irregular  layer. 

The  effects  of  error  in  the  background  bottom  model  upon  the  backscattering 
strength  predicted  by  the  model  is  minimal  because  the  propagation  effects  associated 
with  the  background  bottom  model  is  cancelled  out  in  estimating  the  backscattering 
strength. 

For  the  Gaussian  distribution,  an  acceptable  fit  for  the  upper  irregular  layer  is 
obtained  with  =  ly  =  Im  ,  =  0.2m,  a  —  1%,  and  /?  =  2.  To  us,  these  parameters 

are  less  reasonable.  The  agreements  are  not  as  good  as  with  a  power  law  distribution, 
either.  In  addition,  the  fit  is  far  less  satisfactory  for  the  lower  layer.  All  of  these 
factors  make  the  idea  of  using  a  Gaussian  distribution  to  describe  sound  speed  and 
density  fluctuations  unfavorable  to  us  in  the  ARSRP  backscattering  scenario. 

As  discussed  in  Chapter  3,  for  backscattering  from  random  inhomogeneities,  the 
angular  dependence  of  the  backscattering  strength  is  mainly  determined  by  the  roll¬ 
off  rate  of  its  power  spectrum  versus  wavenumber.  This  is  because  the  roll-off  rate 
is  directly  related  to  the  level  of  small-scale  roughness  at  high  wavenumbers.  For  a 
power  law  distribution,  the  quantity  v  will  control  the  roll-off  rate  {2{i/  +  3/2)  for 
a  3-D  case)  and  therefore  the  angular  dependence.  Different  correlation  lengths  will 
only  alter  the  power  spectrum  level  at  a  particular  wavenumber  but  not  the  roll-off 
rate.  The  aspect  ratio,  defined  by  the  horizontal  correlation  length  over  the  vertical 
correlation  length,  would  play  a  role  as  well,  although  not  as  significant  as  z/,  which 
is  shown  in  Fig.  6-15.  Figure  6-16  shows  the  results  of  the  parameter  study.  The 
two  curves  have  the  same  u,  a,  I3,  same  aspect  ratio  but  different  correlation  lengths. 
It  can  be  seen  that  the  trends  of  the  two  curves  are  very  similar  but  not  exactly  the 
same.  The  differences  can  be  seen  at  450Hz  and  5bOHz,  for  example.  Although  we 
find  the  dash  line  the  best  fit  with  data,  it  is  a  close  call  since  the  differences  are  small 
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Figure  6-15:  Parameter  studies  for  sound  speed  and  density  fluctuations  described  by 
a  power  law  type  of  power  spectrum:  the  effect  of  the  aspect  ratio. 
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Figure  6-16:  Parameter  studies  for  sound  speed  and  density  fluctuations  described 
by  a  power  law  type  of  power  spectrum:  same  aspect  ratio  but  difiPerent  correlation 
lengths,  V  =  0.5. 
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.  This  insensitivity  would  make  it  difficult  to  invert  for  sediment  properties  such  as 
correlation  lengths  from  the  angular  dependence  of  the  backscattering  strength.  On 
the  other  hand,  it  reveals  the  importance  of  estimating  i/  if  random  inhomogeneities 
are  described  by  a  power  law  distribution.  In  the  meantime,  it  would  benefit  the  task 
of  predicting  the  directivity  of  the  backscattered  field  without  full  knowledge  of  the 
ground  truth. 

In  summary,  the  model  is  able  to  match  the  ARSRP  backscattering  data  well  with 
a  power  law  type  of  power  spectrum  description  of  random  sound  speed  and  density 
fluctuations.  Yet,  model/data  comparisons  for  more  data  sets  are  needed  before  any 
conclusion  reached  in  this  section  can  be  generalized.  Nevertheless,  the  successful 
model/data  comparison  does  show  the  merit  of  our  volume  backscattering  model. 
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Chapter  7 


Conclusions  and  Future  Work 


In  this  thesis  an  attempt  has  been  made  to  model  monostatic  backscattering  from 
3-D  volume  inhomogeneities  in  the  sediment  and  to  compare  the  results  of  the  model 
with  ARSRP  backscattering  data.  While  the  work  presented  here  is  helpful  in  under¬ 
standing  the  underlying  physics  and  predicting  the  scattering  phenomena,  more  work 
remains  to  be  carried  out  in  order  to  achieve  one  of  our  ultimate  goals,  solving  the 
inverse  problem.  In  this  chapter,  we  will  try  to  draw  some  general  conclusions  from 
the  work  being  done  and  come  up  with  some  suggestions  for  future  research. 

7.1  Conclusions 

The  results  of  the  thesis  can  be  summarized  as  follows: 

•  A  scattering  process  cannot  be  modeled  correctly  without  properly  accounting 
for  the  incident  field.  The  scattered  field  and  the  incident  field  are  interre¬ 
lated  and  actually  it  is  the  total  field,  i.e.,  the  combination  of  the  two,  that 
controls  the  reradiation  of  acoustic  energy  in  an  inhomogeneous  medium.  In 
some  scenarios,  the  scattered  field  is  much  weaker  than  the  incident  field,  and 
therefore  the  total  field  can  be  approximated  by  the  incident  field  {i.e.,  the 
Born  approximation).  In  any  case,  an  accurate  evaluation  of  the  incident  field 

is  essential  in  modeling  the  scattering  process.  As  revealed  in  Chapter  2,  several 

\ 

approximate  propagation  models  used  in  volume  scattering  studies  sometimes 
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have  their  drawbacks  in  dealing  with  scattering  from  a  bottom  with  complicated 
sound  speed  structure  where  multi-path  effects  are  significant.  As  a  matter  of 
fact,  not  only  do  we  need  to  examine  the  accuracy  of  these  approximation  meth¬ 
ods  in  different  scenarios,  but  also  it  is  often  necessary  to  evaluate  the  incident 
field  exactly  by  either  an  analytical  solution  or  a  numerical  method  such  as 
wavenumber  integration.  The  above  results  are  obtained  through  simulations 
in  a  deep-water  environment.  In  shallow  water,  we  would  expect  the  propa¬ 
gation  of  sound  waves  become  more  intricate  because  of  the  waveguide  effect 
[7],  which  would  further  the  requirement  of  a  good  propagation  model  so  as  to 
calculate  the  incident  field  accurately  in  scattering  modeling. 

•  With  OASES  chosen  as  the  propagation  model  in  this  work,  we  are  provided 
with  the  ability  to  model  volume  scattering  from  a  bottom  with  complicated 
sound  speed  structures  such  as  sound  speed  gradient,  multi-layering,  and  even 
an  elastic  basement  as  long  as  the  scatterers  are  located  in  a  fluid  layer. 

•  The  concept  of  applying  the  equivalent  surface  backscattering  strength  in  char¬ 
acterizing  volume  scattering  processes  is  still  valid  in  some  high-frequency  and 
high-attenuation  bottom  scattering  studies.  However,  precautions  need  to  be 
taken  when  multi-path  contributions  to  the  field  are  evident,  which  would 
put  the  obtained  angular  dependence  of  the  equivalent  surface  backscattering 
strength  in  question. 

•  The  volume  scattering  model  described  in  Chapter  3  provides  a  tool  to  model 
scattering  contributions  from  both  sound  speed  and  density  fluctuations  in  a 
volume.  Simulations  in  Chapter  5  have  confirmed  the  fact  that  density  varia¬ 
tions  are  dominant  in  backscattering.  Meanwhile,  the  analytic  solutions  in  a 
free-space  scenario  show  clearly  that  the  nature  of  the  scattered  field  depends 
strongly  on  the  statistical  distribution  of  the  random  sound  speed/density  vari¬ 
ations.  Both  the  characteristic  length  scale  and  the  power  spectrum  descrip¬ 
tions  of  the  random  inhomogeneities  have  great  impact  on  the  directionality 
of  the  scattered  field.  Although  the  monostatic  configuration  cannot  resolve 
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the  anisotropic  nature  of  the  scattered  field  in  horizontal  planes,  the  angular 
dependence  would  behave  very  differently  because  of  that.  Effects  due  to  the 
vertical  correlation  of  the  random  scatterers  are  evident  in  the  numerical  ex¬ 
periments.  To  model  the  scattering  volume  as  vertically  uncorrelated  stacks  of 
layers  would  bear  the  risk  of  losing  important  contributions  due  to  the  vertical 
correlation.  All  the  above  helps  us  to  reach  the  conclusion  that  it  is  critical 
to  obtain  ground  truth  for  the  sound  speed/density  variations  in  the  sediment 
volume.  The  effort  of  achieving  better  understanding  and  parameterizing  the 
variations  of  the  sediment  properties  should  be  an  integrated  part  of  volume 
scattering  studies. 

•  Simulating  scattering  from  full  3-D  volumetric  inhomogeneities  is  still  very  com- 
putionally  expensive.  In  Chapter  4,  the  monostatic  backscattering  configuration 
has  been  taken  advantage  of  so  that  the  azimuthally-summed  3-D  random  field 
can  be  generated  with  much  relaxed  computational  requirements.  It  enables  us 
to  get  a  grip  on  the  full  scope  of  monostatic  volume  backscattering  in  a  real 
environment. 

•  In  the  model  and  ARSRP  backscattering  data  comparison,  a  power  law  distri¬ 
bution  of  random  sound  speed/density  variations  is  found  to  fit  the  data  very 
well.  Parameter  wise,  on  the  one  hand,  the  horizontal  correlation  length  is  much 
larger  than  the  vertical  correlation  length,  which  is  consistent  with  some  of  the 
geophysical  models  of  the  sediment  (so-called  “pancake”  models).  On  the  other 
hand,  the  change  of  horizontal  correlation  length  in  a  power  law  distribution 
does  not  have  strong  effect  on  the  slope  of  the  backscattering  strength  versus 
grazing  angle  but  only  the  level.  Rather  the  slope  is  determined  by  the  decaying 
rate  of  the  corresponding  power  spectrum  versus  wavenumber,  i.e.,  the  fractal 
dimension,  which  is  directly  related  to  the  sound  speed/density  fluctuations  in 
small  scales.  This  would  benefit  people  who  want  to  predict  scattering  phenom¬ 
ena  but  only  with  limited  knowledge  of  the  sediment  properties.  However,  this 
insensitivity  presents  a  challenge  for  people  who  attempt  to  do  inversion.  All  in 
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all,  this  first-order  volume  scattering  model  has  been  proved  capable  of  match¬ 
ing  the  experimental  data  and  the  power  law  distribution  is  found  to  describe 
the  sound  speed  and  density  fluctuations  in  the  sediment  satisfactorily. 


7.2  Future  Work 

There  are  a  number  of  ways  to  improve  or  extend  the  proposed  model  in  this  thesis. 

•  A  logical  next  step  would  be  to  take  advantage  of  the  ability  rendered  by  OASES 
to  include  rough  surface  scattering  in  the  picture.  The  water/sediment  interface 
scattering  was  insignificant  in  the  ARSRP  backscattering  experiment.  However, 
a  more  general  scenario  would  inevitably  involve  both  rough  surface  scattering 
and  volume  scattering.  It  would  be  interesting  to  see  the  combined  effect  and 
to  compare  that  with  surface  scattering  or  volume  scattering  only. 

•  The  volume  scattering  model  that  we  have  here  is  essentially  a  single  scattering 
model,  which  would  certainly  be  inaccurate  when  multiple  scattering  is  of  con¬ 
cern.  The  development  of  a  multiple  scattering  model  will  improve  our  modeling 
ability  and  further  our  understanding  of  the  scattering  mechanisms. 

•  A  very  useful  extension  of  this  work  would  be  to  model  bistatic  scattering 
instead  of  only  monostatic  backscattering,  although  we  may  not  be  able  to 
enjoy  the  benefit  of  a  simplified  3-D  random  field  generation. 

•  The  work  present  above  models  volume  scattering  in  fluid  media  only  and  cannot 
account  for  scattering  from  an  elastic  bottom.  The  model  would  be  more  useful 
if  it  could  describe  volume  scattering  in  elastic  media,  which  has  attracted  more 
attention  recently. 

•  If  an  experiment  similar  to  the  ARSRP  experiment  were  to  be  performed  in 
the  future,  more  effort  in  gathering  ground-truth  information  is  highly  recom¬ 
mended,  especially  the  attenuation  coeflftcient. 
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Appendix  A 


The  Variance  and  Power  Spectral 
Density  of  the  2-D 
Azimuthally-Summed  Random 
Field  with  Gaussian  Correlation 
Function 


As  mentioned  in  Chapter  4,  the  2-D  Gaussian  correlation  function  in  Cartesian  coor¬ 
dinates  is 

Cfe.y  =a^exp(-i[|  +  |]).  (A.1) 

where  o  is  the  standard  deviation,  and  4  and  ly  are  the  correlation  lengths  in  the  x 
and  y  directions,  respectively.  The  corresponding  power  spectral  density  is 

i-2/2  I  ;l272 

WiK,ky)  =  27rlJyexp{-^^Y'')-  (A.2) 
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For  an  isotropic  random  field,  we  have  lx  —  ly  =  1.  After  the  coordinate  transfor¬ 
mation 

{X  =  rcos<}) 
y  —  rsin<f> 

the  correlation  function  of  e  becomes 


<  e(n,0i)e*(r2,<^2)  >=  exp( 


rj  +  rj-  2rir2Cos((^i  -  (/>2) 
2/2 


), 


(A.4) 


where  (*)  stands  for  the  conjugate  operation.  Since 


r]{r)  =  /  e(r, 
Jo 


the  correlation  function  of  rj  becomes 


exp( 


<  e(ri,(^i)€*(r2,<A2)  > 

I  - p - 


4?r’(T'/o(^)  exp( 


d±il) 

2/2  ^ 


(A.5) 


If  ri  =  r2  =  r,  the  variance  at  range  r  is 


<  v{r)v*{'r)  >=  47rV/o(^)exp(-^). 


(A.6) 


For  an  anisotropic  random  field,  the  correlation  function  of  77  (r)  is 


<  'niri)v*ir2)  > 


r2n  r2n 

/  /  exp(- 

Jq  Jo 


{r2C0S(j)2  —  riCos(f>iy  {r2sin(p2  —  rism^i) 
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+ 


2'^ 


)d<pid^2-  (A.7) 
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The  variance  as  a  function  of  r  in  this  case  becomes 


Var{'q{r)}  = 


^  t2n  r27r  ^2 


The  power  spectral  density  is  as  follows 


W,(K)  = 


=  2nlJ,exp{-^^!^^)£’ 

=  4.^y,exp(-^,l(M), 


(A.9) 
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Appendix  B 


Estimation  of  the  Standard 
Deviation 


The  standard  deviation  is  widely  used  to  describe  the  dispersion  of  a  random  distribu¬ 
tion.  It  is  known  to  work  best  for  a  Gaussian  distribution.  However,  the  values  of  the 
standard  deviation  is  greatly  affected  by  outliers  for  some  non-Gaussian  distributions 
and  therefore  may  be  a  poor  measure  of  dispersion  in  those  cases.  The  formula  to 
estimate  the  standard  deviation  from  random  realizations  is 


STD{y}  = 


^|  AT-l 


N 


-  2/)^ 


where  N  represents  the  number  of  realizations  and  y  the  mean  of  y. 


(B.l) 


For  the  simulated  backscattering  coefficients  from  the  model  in  Chapter  6  (a  power 
law  distribution),  the  distribution  is  not  Gaussian  but  more  like  an  exponential  dis¬ 
tribution  (or  a  Chi-square  distribution  with  2  degrees  of  freedom),  which  can  be  seen 
in  Fig.  B-1.  It  is  clear  that  in  this  case  the  standard  deviation  no  longer  qualifies  as 
a  quantitative  measurement  of  the  dispersion. 


The  backscattering  coefficient  estimated  from  the  data  is  the  mean  of  eight  pings 
(eight  independent  realizations).  If  the  mean  value  is  treated  as  a  new  random  vari- 
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Figure  B-1:  Histogram  of  the  simulated  backscattering  coefficients  at  0  =  42  degrees 
and  /  =  500i/z  (200  realizations). 
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Backscattering  coefficient  averaged  over  8  realizations  x  10"^ 


Figure  B-2:  Histogram  of  the  simulated  backscattering  coefficients  averaged  over  8 
realizations  at  0  =  42  degrees  and  /  =  hQQHz  (200  realizations). 

able,  we  will  arrive  at  a  different  random  distribution.  Assuming  that  the  backscat¬ 
tering  coefficient  before  averaging  possesses  an  exponential  distribution,  the  new  ran¬ 
dom  variable  will  have  a  Chi-square  distribution  with  16  degrees  of  freedom.  In  other 
words,  it  is  closer  to  a  Gaussian  distribution  than  the  original  exponential  distribu- 
tion(it  will  be  exactly  a  Gaussian  distribution  if  the  average  is  performed  over  an 
infinite  number  of  independent  realizations  according  to  the  Central  Limit  Theorem). 
This  can  be  confirmed  by  the  histogram  of  the  simulated  backscattering  coefficient 
averaged  over  eight  realizations  as  shown  in  Fig.  B-2.  As  a  result,  the  standard 
deviation  estimated  from  the  above  equation  would  be  a  much  better  measurement 
of  the  dispersion,  which  is  shown  as  the  error  bars  in  Fig.  6-13  and  Fig.  6-14. 

While  for  the  data,  only  8  pings  are  available.  Therefore,  the  above  procedure 
done  to  200  realization  of  simulated  backscattering  coefficient  cannot  be  performed, 
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o> 


Figure  B-3:  The  estimated  backscattering  strength  with  standard  deviations  (error 
bars)  over  all  the  frequencies  for  the  (a)  upper  and  (b)  lower  irregular  layers  at  the 
east  side  of  the  sediment  pond. 

which  is  the  reason  why  error  bars  are  not  plotted  in  Figs.  6-13  and  6-14  to  rep¬ 
resent  the  standard  deviation  estimated  from  the  data.  However,  the  average  over 
independent  frequencies  in  each  ping  should  provide  results  which  are  similar  to  those 
obtained  by  averaging  over  pings.  By  using  this  approach,  we  can  obtain  reasonably 
good  estimates  of  the  standard  deviation  of  the  backscattering  strength  from  the  data, 
which  are  shown  in  Fig.  B-3. 
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